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Abstract 

A  procedure  is  described  for  developing  simple 
approximate  equations  of  state  of  liquids  from  Hugoniot  P-V 
relations  determined  in  shock  wave  measurements.  This  is 
applied  to  a  number  of  liquids  and  a  table  of  coefficients  is 
given. 

The  formalism  of  irreversible  thermodynamics  is 
applied  to  time-dependent  phase  transitions  in  iron  and  an 
approximate  set  of  constitutive  relations  is  obtained  in  a 
form  suitable  for  numerical  integration  with  the  equations  of 
continuum  dynamics.  These  are  applied  in  an  approximate  form 
to  study  the  development  of  the  two-wave  structure  in  iron 
caused  by  the  or-e  phase  transition. 

Finite  strain  theory  is  applied  to  the  analysis  of 
shock  wave  data  for  quartz,  and  the  results  supply  enough 
information  to  estimate  some  of  the  fourth-order  elastic 
constants. 
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Foreword 


The  work  reported  here  represents  the  results  and  the 
state  of  understanding  which  had  been  achieved  in  December, 

1966.  Since  that  time  some  further  progress  has  been  made, 
and  it  will  be  described  in  later  reports  of  this  series. 
Particular  attention  is  drawn  to  the  approximate  temperature 
calculation  described  in  Section  4.1.  The  basis  of  this  approx¬ 
imation  is  now  thought  to  be  unsound,  so  some  skepticism  should 
be  maintained  concerning  temperature  effects  reported.  For¬ 
tunately  these  are  few  and  slight,  and  the  general  conclusions 
of  the  report  are  not  affected  by  possible  errors  here. 

A  major  shortcoming  of  the  theory,  but  one  which  is  very 
hard  to  evaluate,  is  the  complete  reliance  on  equilibrium 
thermodynamics  to  describe  the  static  behavior.  It  is  quite 
likely  that  metallurgical  considerations  govern  the  true  pro¬ 
gress  of  the  <*-€  transition  in  Iron;  consequently  the  static 
reference  states  to  which  dynamic  effects  are  referred  may  be 
metastable  and  very  different  from  the  thermodynamic  states. 

This  is  suggested  by  recent  static  pressure  measurements  by 
Bassett  and  co-workers  (J.  Appl.  Phys.,  Jan.  1967)  and  by  shock 
de-magnetization  experiments  by  R.  A.  Graham  of  Sandia  Corpora¬ 
tion  (private  communication) .  Satisfactory  metallurgical 
models  for  such  processes  are  not  presently  available;  it  is 
hoped  that  shock  wave  measurements  will  help  to  stimulate  the 
development  of  such  models. 

February,  1968 
xi 


I.  INTRODUCTION 


Numerical  methods  for  the  integration  the  equations 
of  finite  amplitude  wave  propagation  have  reached  such  a  stage 
of  development  that  the  principal  limitation  on  predictions  of 
wave  effects  is  more  apt  to  be  uncertainty  in  the  constitutive 
relations  of  the  material  than  inability  to  perform  the  inte¬ 
grations.  The  complete  constitutive  relations  for  a  real  solid 
may  at  present  be  regarded  as  unknowable.  Practically  useful 
relations  can  be  obtained  from  a  succession  of  approximations, 
to  each  of  which  is  attached  some  uncertainty.  For  high  ampli¬ 
tude  compressive  waves  the  predominant  relation  is  the  hydro¬ 
static  one  between  pressure,  volume  and  temperature.  To  this 
may  be  added  the  effects  of  finite  shear  dtrength,  which  makes 
the  pressure  tensor  anisotropic,  strain-rate  effects  which  cause 
deviations  from  equilibrium,  and  partition  of  internal  energy 
among  thermal,  surface,  and  inhomogeneous  effects. 

In  the  present  work  the  primary  emphasis  is  on  the 
effects  of  phase  transitions  on  compressive  wave  forms,  partic¬ 
ularly  when  the  rate  of  transition  is  too  slow  to  exactly  follow 
the  changes  in  pressure,  temperature,  and  density  associated 
with  the  compressive  wave  which  initiates  the  transition.  In 
order  to  study  these  effects  it  has  been  necessary  to  develop 
computer  programs  for  integration  of  the  flow  equations  for  the 
appropriate  constitutive  relations.  In  the  course  of  this 
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development  other  useful  results  on  equations  of  state  have  been 
produced,  and  these  are  described  in  the  following  sections. 


II.  CONSTITUTIVE  RELATIONS 


2.1  General  Considerations 

The  term  "constitutive  relations"  is  used  in  a  generic 
sense  to  encompass  all  material  properties  which  must  be  combined 
with  the  equations  of  continuity,  motion,  and  energy  conserva¬ 
tion  to  supply  a  complete  set  of  flow  equations.  Constitutive 
relations  in  practice  usually  reduce  to  an  equation  of  state 
relating  pressure,  volume,  and  temperature  or  internal  energy. 
This  simplification  is  partially  enforced  by  ignorance  of  other 
material  relations;  it  often  yields,  in  addition,  quite  a  good 
description  of  wave  propagation  over  a  wide  range  of  parameters. 
The  eauation  of  state  is  necessarily  accompenied  by  a  statement 
about  the  variation  of  specific  heat  with  pressure  and  tempera¬ 
ture.  It  may  on  occasion  include  information  about  rigidity  and 
yield  or  even  rate  effects.  For  the  present  we  consider  the 
eauation  of  state  as  defined  above.  These  considerations  are 
themselves  useful  and  they  will  provide  insight  into  the  require¬ 
ments  which  must  be  satisfied  by  more  general  constitutive  rela¬ 
tions.  The  following  remarks  are  necessarily  limited  in  their 
scope.  More  detailed  information  will  be  found  in  various  review 
articles  (1,2). 

2.2  Equations  of  State  for  Fluids  (3) 

A  "complete"  equation  of  state  for  a  fluid  is  a  rela¬ 
tion  between  thermodynamic  variables  which  is  sufficient  for 
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calculating  any  thermodynamic  parameter  of  the  material,  given 
two.  For  example,  if  specific  internal  energy,  E,  is  given  as 
a  function  of  specific  entropy,  S,  and  specific  volume  v, 

E  =  E(S,v) ,  then 

p  =  -(*e/*v)s 

T  =  (*e/*S)v 

H  =  E  +  pv 
etc . , 

where  p,  T,  H  are  pressure,  temperature  and  specific  enthalpy. 
On  the  other  hand,  if  pressure  is  given  as  a  function  of  T  and 
v,  as  usually  occurs,  then  neither  internal  energy  nor  enthalpy 
can  be  calculated  without  specifying  the  specific  heat. 

When  partial  equations  of  state  are  given,  as  in  the 
last  example,  then  certain  limitations  are  placed  on  other 
thermodynamic  quantities  if  all  are  to  be  compatible.  This  is 
illustrated  by  the  following  example.  Suppose  that  specific 
heat  at  constant  pressure  is  known  as  a  function  of  temperature 
and  assumed  to  be  a  function  of  temperature  alone:  Cp  =  Cp(T) . 
Then  by  the  following  argument  we  can  see  that  the  relation  be¬ 
tween  p,  v,  and  T  must  be  that  of  Eq.  (2.1)  belox^: 

(aH/*T)p  =  Cp(T)  by  definition.  (2.1) 

If  H  »  H(p,T) ,  then 

dH  =  CpdT  +  (*H/*p)T  dp 
or 

H  =  J  Cp(T)dT  +  f(p)  . 
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We  may  thus  write  that 

(*R/*p)T  =  f'(p)  =  (ftH/*p)s  +  (*HAs)p(*S/*p)T 

=  v  -  T  (*v/*T)p. 

2 

Dividing  by  T  we  find  that 

f’(p)/T2  =  -(*(v/T)/?'T)p  . 

When  integrated  thi~  yields  the  implied  form  for  the  equation 
of  state: 

v  =  f'(p)  +  Tg(p)  (2.2) 

where  g(p)  =  (W*T)p 

In  a  similar  way  it  can  be  shown  that  if  specific  heat  at  con¬ 
stant  volume,  Cv,  is  a  function  of  T  alone,  Cv  88  Cy  (T) ,  then 
it  must  follow  that  p,  T,  and  v  are  related  by  a  Gruneisen 
equation: 

p  =  f(v)  +  Th(v)  (2.3) 

3oth  Eqs.  (2.2)  and  (2.3)  are  remarkably  simple.  The 
(p,  v,  T)  surfaces  that  can  be  represented  by  either  of  these 
can  be  constructed  from  a  "baraboo-place-mat. "  In  the  first 
case  the  straight  sticks  lie  ir  planes  of  constant  p;  in  the 
second  they  lie  in  planes  of  constant  v. 

Eq.  (2.3)  has  the  same  form  as  the  Mie-Gruneisen  equation: 

P  -  Pk(v)  +  TCV  (T-Tq)/v  (2.4/ 

where  pk  is  pressure  on  the  isotherm  T  =  T0.  With  the  assumption 
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that  Cv  =  CV(T),  Eqs.  (2.2)  and  (2.3)  then  imply  that 

(*(rcv)/*T)v  =  0. 

This  condition  implies  that  (ar/*T)v  0  unless  Cv  =  const. 

While  it  is  very  likely  true  that  r  does  depend  upon  T  (4) , 
the  only  available  theories  of  any  generality  suppose  that 
T  =  F(v)  (1).  Fowles  has  obtained  a  more  general  compatibility 
relation  for  r  and  Cv  (5) : 

(*Cv/*lnv)T  =  (^(rCv)/ftlnT)v 

It  turns  out  that  this  is  satisfied  by  the  Debye  theory  of 
specific  heats.  However,  any  dependence  of  Debye  temperature 
on  temperature  violates  the  condition. 

Despite  these  difficulties,  Eq.  (2.3)  has  been  commonly 
used  to  extend  pressure-volume  data  determined  from  shock 
studies  into  off-Hugoniot  regions  (1).  Doran  (6)  has  gone  even 
farther  to  show  that  quite  reasonable  representations  of  the 
equations  of  state  of  solids  can  be  obtained  with  r/v  =  constant. 

In  view  of  the  limitations  on  r  and  Cv,  it  is  unlikely 
that  the  zero  degree  isotherms  calculated  using  Eq.  (2.4)  and  the 
Slater  or  Dugdale-McDonald  relation  for  r  are  physically  reli¬ 
able.  Their  principal  virtue  is  that  they  provide  a  consistent 
and  reproducible  procedure  for  determining  a  reference  curve. 

An  alternative  procedure  is  to  calculate  the  isotherm  passing 

through  the  initial  state  of  the  Hugoniot.  This  avoids  some  of 

o 

the  difficulties  associated  with  extrapolating  to  0  K.  It  in¬ 
troduces  some  new  ones  inasmuch  as  there  is  now  no  theory  for 
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calculating  r.  However,  since  the  low  temperature  region  is 
eliminated,  the  assumption  r/v  =  constant  may  be  a  reasonable 
one.  Then  Eq.  (2.4)  becomes 


p  =  Pi  +  b  (E-E^) 

^  =  Eo  +  %  Pr  (Vv> 

dEi  =  (To  C^p/^T)v  -  pt)dv 


(2.5) 

(2.6) 
(2.7) 


=  (bCvT0  -  Pi)dv 

where  b  =  r/v  =  constant,  Cv  =  constant,  p^(v)  and  E^(v)  are 
pressure  and  internal  energy,  respectively,  on  the  Tq  isotherm, 
and  subscript  "h"  refers  to  the  Hugoniot  curve.  Setting  p  and 
E  in  Eq.  (2.5)  equal  to  p^  and  E^  and  combining  with  Eqs.  (2.6) 
and  (2.7)  yields  a  differential  equation  for  p^: 

(dp^/ dv)  +  bp^  =»  (1  -  b(vQ~v)/2)(dpH/dv) 

+  bpH/2  +  b2CvTQ  (2.8) 

The  solution  of  this  equation  is 


P^v)  =  A  exp  (-bv)  +  bCvTQB  (2.9) 

A  =  f(v)  -  bf  f(v)dv 
v  o 

B  =  1  -  exp  (b(vQ-v)) 


f(v) 


(1  “  (b/2)(v0-v))  pH  exp  (bv). 
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Experience  has  shown  that  Hugoniot  data  for  liquids  and  solids 
can  be  fitted  quite  well  by  curves  of  the  form 

3 

Ph(v)  “  I  anxT1  (2.10) 

n=l 

where  x  =*  Pv  -  1. 

o 

Equations  of  the  form  (2.10)  have  been  fitted  to  shock 
data  on  liquids  and  used  to  calculate  p^(v)  and  E^(v)  from  Eqs. 
(2.9)  and  (2.7),  respectively.  The  numerical  results  are  used 
in  a  least  squares  procedure  to  calculate  the  coefficients  bfl 
in  the  equation  for  isothermal  pressure: 

Pi  ■  I  vn 

n=l 

where  x  =  Pv0“l  as  in  Eq.  (2.10).  The  coefficients  afl  and  bR 
are  given  in  Table  I. 
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2.3  Elastic  Solids  (G.  R.  Fowles) 

The  shock  compression  of  quartz  is  of  particular  interest 
because  of  its  importance  to  geophysics,  its  wide-spread  use  in 
shock  wave  studies  as  a  pressure  transducer,  and  because  it  rep¬ 
resents  a  different  class  of  materials  from  the  more  thoroughly 
studied  metals.  In  this  paper  we  describe  measurements  similar 
to  those  reported  by  ijackerle  (15").  The  dcta  are  in  substantial 
agreement;  however,  the  recording  techniques  were  somewhat  dif¬ 
ferent  so  that  the  present  results*  provide  independent  corrob¬ 
oration,  in  most  respects,  of  Uackerle's  data. 

In  addition  to  describing  the  experiments  and  the  results, 
we  examine  the  agreement  between  the  uniaxial  stress-strain  data 
derived  from  shock  experiments  and  predictions  based  on  finite 
strain  theory  and  the  second  and  third-order  elastic  constants 
measured  by  McSkimin,  et  al.  (39),  and  Thurston,  et  al.  ( 40)  • 

From  this  comparison  it  is  clear  that  shock-wave  measurements  and 
low  pressure  acoustic  measurements  are  complementary  methods  for 
evaluating  higher  order  elastic  coefficients. 

In  Section  2.31  we  describe  the  experimental  technique 
and  the  experimental  results;  Section  2.32  gives  a  brief  outline 
of  finite  strain  theory  and  its  application  to  the  shock  experi¬ 
ments.  Conclusions  are  discussed  in  Section  2.33. 


*These  data  were  reported  originally  in  the  author's  Ph.D. 
thesis  (48) . 


2.31  Experiments 
A.  Experimental  Method 
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In  the  experiments  shock  propagation  velocities  and  associated  free 
surface  velocities  were  measured  in  alpha  quartz  crystals  oriented  as  X,  Y, 
or  Z-cuts.*  Shock  waves  of  varying  intensity  were  generated  by  plane-wave 
explosive  lenses  with  or  without  additional  explosive  pads. 

The  experimental  arrangement  is  shown  schematically  in  Fig.  2.1.  A 
four-inch  diameter  explosive  lens  (and  in  some  cases  an  explosive  pad)  was 
cemented  to  one  surface  of  a  1 /2-inch  thick,  5-inch  diameter  Dural  plate.  The 
quartz  specimens  (usually  two)  were  cemented  to  the  opposite,  lapped  surface 
of  the  plate.  The  specimens  were  accurately  flat  and  polished;  the  tolerance 
on  crystallographic  orientation  was  *  1°.  The  faces  of  the  specimens  in 
contact  with  the  plate  were  vapor-plated  with  aluminum  to  yield  a  reflecting 
surface.  Lucite  mirrors,  also  aluminized  on  their  inside  faces  were  cemented 
to  the  outer  surfaces  of  the  specimens  at  angles  of  3  to  8°.  The  edge  of  the 
lucite  mirror  in  contact  with  the  specimen  was,  in  each  case,  set  back  from 
the  edge  of  the  specimen  at  least  one  specimen  thickness  to  avoid  interference 
from  edge  effects.  The  X,  Y,  or  Z  orientation  refers  to  the  smallest  linear 
dimension  and  also  designates  the  direction  of  shock  propagation.  The  X-cut 
crystals  were  measured  in  both  the  +  and  -  orientations  because  of  the  large 
differences  observed  in  electrical  experiments  (41). 

In  some  of  the  experiments  an  inclined  lucite  mirror  was  cemented 
directly  to  the  aluminum  plate.  Its  function  was  to  measure  the  free-surface 
velocity  of  the  aluminum  to  permit  impedance-match  solutions  to  the  final 
shocked  statas  (1). 


*Synthetic  crystals  supplied  by  Valpey  Corporation. 
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ALUM  I  NUM 
PLATE 

CAMERA  SLIT 


(VACUUM  CHAMBER 
NOT  SHOWN) 


TOP  VIEW 


Fig. 2.1  Diagram  of  experimental  assembly 
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Tna  angles  of  the  inclined  mirrors  with  respect  to  the  quartz  surfaces 
ware  measured  after  assembly  by  mounting  the  assembly  on  a  mill  table  and 
observing  with  a  telescope  the  superposition  of  a  cross-hair  and  its  image 
reflected  alternately  in  the  quartz  and  lucite  surfaces.  The  angles  could 
tr.u  be  measured  to  a  precision  of  0.1%.  Soma  cifficulty  was  encountered  in 
keeping  t:.e  lucite  mirrors  extremely  flat.  It  was  necessary  to  allow  angular 
deviations  of  up  to  *  one  minute  of  arc.  In  each  case  this  amounted  to  less 
than  1/2%  of  the  total  angle. 

In  order  to  obtain  the  desired  accuracy  in  shock  velocity,  1  1%, 
good  contact  (0.0002  inch)  between  the  inside  edge  of  the  inclined  mirror  and 
toe  cuter  quartz  surface  was  required.  A  contact  such  that  r.o  transmitted 
light  was  visible  was  considered  satisfactory. 

In  order  to  avoid  camp! i cations  due  to  air  shocks  the  assembly  was 
evacuated  prior  to  firing  to  a  pre^^ure  of  less  than  0.05  torr.  A  hemicylindr 
section  of  'ucite  tubing  cemented  to  the  aluminum  plate  served  as  a  vacuum 
chamber. 

A  pnotograph  of  an  assembly,  without  explosive,  prior  to  firing  is 
shown  as  Fig.  2.2. 

The  assembly  was  viewed  through  a  slit  of  a  rotating  mirror  streak 
camera  aligned  along  the  centers  of  the  inclined  mirrors  in  the  direction  of 
maximum  inclination  (i.e.,  the  direction  in  which  the  mirror  angles  were 
previously  measured).  The  slit  width  was  0.05  mm;  the  time  resolution, 
determined  from  the  slit  width  and  the  camera  writing  speed  (*3.h1  mm/ys),  was 
approximately  0.01  ys. 

II -.jmi nation  was  provided  by  an  explosive  argon  light  source  consisting 
of  a  4-incn  diameter,  18-inch  long  cardboard  tube  with  a  one* inch  pad  of 
composition  C-3  explosive  at  one  end.  A  ground  glass  diffusing  screen  was 
placed  over  the  other  end  and  argon  was  flowed  through  the  tube  continuously. 
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Fig.  2.2. — Photograph  of  Experimental  Assembly 


I 


]  5 

The  light  source  explosive  was  initiated  simultaneously  with  the  plane  wave 
lens  of  the  experimental  assembly;  the  resulting  strongly  luminuous  shock  in 
the  argon  produced  a  bright  reflection  from  the  aluminized  surfaces  a  few 
microseconds  before  the  first  arrival  to  be  recorded  in  the  quartz. 

A  drawing  of  the  complete  arrangement  as  it  appeared  before  firing  is 
shown  as  Fig.  2.3. 

An  abrupt  change  in  intensity  of  the  light  reflected  from  the  aluminized 
surfaces  of  the  assembly  showed  arrival  times  of  the  shock  fronts  and  free 
surfaces  upon  impact  with  the  mirrors. 

A  streak  camera  photograph  taken  in  this  manner  is  shown  in  Fig.  2.4. 

The  two  specimens  in  this  shot  were  Z-cut;  the  upper  one  was  1/8-inch  thick 
and  the  lower  1/4-inch  thick.  The  final  pressure  was  approximately  200  kbar. 

Az  time,  Tg,  the  reflection  from  the  rear  (aluminized)  face  of  the  quartz 
extinguishes  abruptly  as  the  shock  arrives  at  the  quartz-aluminum  interface. 

At  time,  T-j,  the  first  shock  arrives  at  the  quartz  free  surface.  The  traces 
are  relatively  smooth  until  the  change  in  slope  caused  by  the  arrival  of  the 
second  shock  at  time  T2;  thereafter  the  traces  are  slightly  irregular.  A 
slight  curvature  to  the  trace  of  the  first  shock  can  be  detected.  This  slowing 
up  of  the  free-surface  is  due  to  stress-relaxation  effects,  as  was  pointed  out 
by  Wackerle(15)  , 

For  reliable  results  the  point  of  collision  of  the  quartz  free  surface 
with  the  inside  surface  of  the  mirror  must  travel  with  supersonic  velocity 
with  respect  to  both  quartz  and  lucite  (non-jetting  configuration).  Consequently, 
the  initial  mirror  angle  must  be  less  than  approximately 

a  max  =  sin"^  uf 

where  u^  is  the  quartz  free-surface  velocity  and  is  the  larger  of  the  two 
shock  wave  velocities  in  quartz  and  lucite.  This  criterion  restricted  the 

usable  mirror  angles  to  less  than  about  8°. 
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B.  Data  Reduction 

The  shock  velocities  were  determined  from  distances  measured  on  the 
film  and  the  known  writing  speed  of  the  camera.  The  velocity  of  the  second 
shock  requires  corrections  because  of  the  motion  of  the  free  surface  and 
because  of  the  interaction  of  the  second  shock  with  the  reflection  of  the 
first  shock.  The  first  of  these  is  straightforward  and  a  simple  derivation 
gives : 


d  +  ufl(T2  -  T-j) 
u2  "  ' 


(2.11) 


T2  "  T0 

where  d  is  _ne  initial  specimen  thickness,  u^-j  is  the  free-surface  velocity 
due  to  the  first  shock,  and  Tq,  T-j  ,  and  T2  are  the  arrival  times  of  the  shock 
fronts  as  shown  in  Fig.  2.4. 

The  correction  due  to  the  interaction  of  the  second  shock  with  the 
reflection  of  the  first  requires  knowledge  of  the  state  (and  constitutive 
relation)  of  the  quartz  in  the  region  between  the  two  fronts  and  cannot  be 
made  unequivocally.  However,  the  assumption  that  the  material  is  stressed  and 
relieved  only  elastically  by  the  first  wave  leads  to  a  large  correction  and 
unreasonably  high  compression  for  the  state  behind  the  second  shock  in  shot 
No.  7394  (Table  II)  The  results  from  that  shot  are  the  most  sensitive  to  this 
correction  because  the  second  shock  was  relatively  slow  with  respect  to  the 
first.  For  the  other  experiments  the  correction  is  smaller  and  does  not 
appreciably  affect  the  conclusions. 

It  should  be  emphasized,  however,  that  the  result  for  shot  7394  implies 
that  an  irreversible  change  in  the  material  properties  occurs  between  the  two 
shock  fronts.  This  conclusion  is  consistent  with  the  observed  relaxation  of 
the  state  of  the  first  shock.  It  is  not  consistent  with  an  assumption  of 
conventional  elastic-plastic  behavior. 

Because  of  the  arbitrariness  regarding  the  interaction  correction  the  data 
are  here  reported  without  such  a  correction.  The  correction  used  by  WACKERLE  (15) 
is  plausible,  but  does  not  significantly  change  the  data  of  this  paper. 


distance  along  slit 


Fig.  2. 4. --Streak  Camera  Photograph  Showing  Shock 
Arrival  Times  and  Free-Surface  Traces.  Shock  No.  7394. 


SUMMARY  OF  EXPERIMENTAL  DATA 


|  Z  3.411  0.462  0.631  .1.51  2.47  7.33  0.751  147.6  O.oSill  67.7  5.49  1.23  215.3  0.U26  211.2 

[7395  P-40  Y  6.601  1.08S  1.448  0.836  1.60  6.07  0.418  67.7  0.9300  22.2  4.77  0.000  114.9  0.8195  8S.3 

!  Y  3.399  0.562  0.745  0.862  1.58  6.05  0.431  69.3  0.9287  22.2  4.77  0.790  113.9  0.8319  85.5 
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T.'.a  free-surface  velocities  were  calculated  from  the  measured  slopes 


of  the  traces  by  means  of  the  relation: 


u-  =  tan  a1 
T  HF  tan  Y' 


(2.12) 


where  a'  is  the  effective  angle  of  the  inclined  mirror  with  respect  to  the 
eiartz  surface,  y'  is  the  angle  of  the  trace  on  the  film  wian  respect  to  the 
space  ax*.s,  M  is  the  magnification  or  ratio  of  distance  on  the  film  to  the 
corresponding  distance  on  the  shot,  and  F  is  the  writing  speed  of  the  camera. 
The  parameters,  a',  and  y',  of  this  relation  are  not  identical  to  their  nominal 
values,  ^  and  y  because  of  tilt  of  the  incident  shock  and  slight  departures 
from  orthogonality  of  the  slit  and  sweep  directions.  The  corrections  are 


given  by 


tan  a'  =  tan  a  (1  +  e'/tan  y) 


tan  y'  =  tan  y  sec  5  (1  -  tan  y  tan  s) 


where  a  is  the  angle  of  the  inclined  mirror  with  respect  to  the  quar-z  surface, 
6*  is  the  angle  of  shock  tilt  as  measured  on  the  film,  <$  is  the  angle  of  the 
slit  with  respect  to  the  normal  to  the  sweep  direction,  and  y  is  the  angle  of 
the  trace  with  respect  to  the  slit  direction  (Fig.  2.5). 

Tr.e  observed  shock  wave  velocities  and  associated  free-surface 
velocities  are  given  in  Table  H,along  with  the  initial  conditions  for  each 
experiment  and  ucher  quantities  derived  from  the  measured  velocities. 

Tr.e  experimental  precision,  based  on  assembly  tolerances,  camera 
resolution,  and  film  reading  errors  is  estimated  to  be  *1%  in  shock 
velocity  and  ±5£  in  free-surface  velocity.  Most  of  the  error  in  free- 
surface  velocity  is  due  to  uncertainty  in  reading  the  angle  y * ( *  1°). 
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SHOT 


FILM 


GA-PLTR-003-6I-44A 


Fig.  2.5. — Definition  of  Parameters  Used  in  Adjusting 
Streak  Camera  Data. 
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C.  Experimental  Results 

The  observed  shock  velocities  are  plotted  as  functions  of  the  shock 
particle  velocities  (taken  to  be  one  half  the  free-surface  velocities)  in 
Fig. 2. 6.  Data  by  WACKZRLE  (15)  and  6REGS0N  (42)  are  also  shown.  The 
data  of  Wackerle  are  his  "average"  values,  shown  for  comparison  because  they 
were  determined  or  the  same  basis  as  the  present  results.  The  solid  curves 
are  predicted  from  finite  strain  theory,  to  be  discussed  inSection  2. 33. 

The  agreement  among  the  experimental  data  is  seen  to  be  generally 
satisfactory.  The  only  significant  disagreement  occurs  for  Z-cut  crystals 
at  a  particle  velocity  of  1.23  mm/ys.  The  source  of  this  discrepancy  is 
unknown.  A  shot  fired  by  Gregson  to  remeasure  this  state  agrees  better  with 
tr.e  present  data*.  However,  if  the  present  data  are  correct  (rather  than 
Wackerle's)  some  anomalous  behavior  is  noted  in  the  pressure-volume  plane, 
as  discussed  below. 

The  stress -comp res si  on  states  were  calculated  from  the  measured 
velocities  by  means  of  the  Rankine-Hugoniot  jump  conditions  (i): 


V/VQ  =  1 


UI  -  u0 
UI  "  u0 


al  -  °0  55  P0^UI  '  u0)(uI  '  uo) 

In  these  equations,  V  is  specific  volume,  u  is  particle  velocity, 

U  is  shock  velocity,  a  is  stress  normal  to  the  shock  front,  ar.d  p  is 
density.  Subscripts  0  refer  to  the  state  ahead  of  the  shock;  subscripts  I 
refer  to  t\a  state  behind  the  shock.  Velocities  are  with  respect  to 
laboratory  coordinates.  Figs.  2.7,  2.8,  and  2.9  showthe  results  inthe  stress, 
speci fi c-vcl ume  plane  for  X,  Y,  and  Z  crystals  respectively.  Bridgman's 


*Shown  in  fig.  2.9. 


SHOCK  VELOCITY,  Ut,  mm/fJis 
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Fig.  2.6. — Shock  Velocity  as  Function  of  Particle 
Velocity.  Curves  labelled  3rd,  4th,  are  fits  based  on  zero- 
pressure  elastic  constants  up  to  3rd  and  4th  order  respec¬ 
tively  for  X  and  Z-cut  crystals. 


0.72  0.76  0.80  0.84  0.88  0.92  0.96  1.0 

RELATIVE  VOLUME  V/\£ 


Fig.  2.7. — Stress-Volume  States  Resulting  From  Shock 
Compression  of  X-cut  Quartz.  Solid  curve  is  Bridgman's 
Hydrostatic  Data.  Curves  Labelled  3rd,  4th,  are  fits 
based  on  zero-pressure  elastic  constants  to  3rd  and  4th 
order . 


STRESS -a*- KILOBARS 
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Fig.  2. 8. --Stress-Volume  States  Resulting  From  Shock 
Compression  of  Y-cut  Quartz.  Solid  curve  is  Bridgman's 
Hydrostatic  Data. 


STRESS  — cr-KILOBARS 


z-  CUT 

A-FOWLES 
□  —  GREGSON 
A  -  WACKERLE 


BRIDGMAN 


0.72  0.76  0.80  0.84  0.88  0.92  0.96  1.0 

RELATIVE  VOLUME  V/VQ 

Fig.  2.9 .--Stress-Volume  States  Resulting  From  Shock 
Compression  of  Z-cut  Quartz.  Solid  curve  is  Bridgman’s 
Hydrostatic  Data.  Curves  labelled  3rd,  4th,  are  fits  based 
on  zero-pressure  elastic  constants  to  3rd  and  4th  order. 

The  curve  labelled  X  represents  the  tangential  stresses, 
calculated  from  constants  up  to  and  including  3rd  order. 
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hydros  ta  vie  curve  based  on  measurements  to  98  kbar  is  shown  for  comparison 
(43).  The  curves  labelled  3rd,  4th,  are  fits  based  on  low 

pressure  acoustic  measurements  and  finite  elastic  strain  theory  (Section  2.33). 
These  plots  show  clearly  the  important  features  of  the  compression, 

namely: 

(1)  extremely  high  amplitude  elastic  waves,  up  to  150  kbar 
in  Z-cut  quartz 

(2)  loss  of  rigidity  above  the  elastic  limit,  as  shown  by 
the  agreement  of  the  higher  pressure  shock  data  with 
extrapolation  of  the  hydrostatic  data 

(2)  lack  of  a  unique  value  for  the  Hugoniot  elastic  limit. 

This  behavior  implies  that  yielding  is  not  due  to  dislocation  motion 
as  in  a  metal,  but  is  analogous  (or  identical)  to  fracture.  It  is  shown 
below  that  the  shear  stresses  behind  the  elastic  shocks  approach  the  theo¬ 
retical  snear  strength  of  the  crystal  lattice. 

Tr.e  range  of  the  present  data  is  not  sufficient  to  show  clearly  the 
transformation  to  stishovite,  as  indicated  by  Wackerle's  higher  pressure 
measurements. 

2.33  Finite  Strain  Theory 

Because  the  strains  behind  the  elastic  shocks  are  relatively  large, 
it  is  of  interest  to  examine  the  agreement  of  the  data  with  predictions  of 
finite  strain  theory.  Predictions  are  made  possible  by  the  work  of 
THURSTON  (40)  and  McSKIMIN  (39)  and  their  co-workers  on  the  third-order 
elastic  constants  of  quartz.  Such  comparisons  should  indicate  the  extent 
to  which  third-order  constants  are  sufficient  to  describe  the  stress-strain 
behavior  at  strains  of  the  order  cf  5  -  10%.  The  constants  are  determined 
from  precise  acoustic  measurements  at  strains  of  less  than  0.1%.  ANDERSON  (44) 
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has  already  shown  that  the  second  and  third-order  constants  alone  give 
reasonably  good  predictions  for  hydrostatic  compressions  of  up  to  about 
!5ii  in  quaraz,  provided  the  constants  are  used  in  the  Murnaghar.  logarithmic 
equation  or  the  Birch  equation  of  state. 

Discrepancies  between  the  observed  and  predicted  stress -strain 
curves  can  ce  used  alternatively  to  evaluate  fourth  and  higher  order  con¬ 
stants,  or  to  examine  the  effects  of  adopting  alternate  definitions  of 
strain,  as  suggested  by  KNOPOFF  (45).  Finally,  to  the  extent  that  the 
third-order  constants  give  adequate  predictions,  the  stresses  tangential  to 
the  shock  fronts  can  be  calculated  from  the  observed  stresses  r.ormal  to  the 
fronts  and,  hence,  the  shear  stresses  sustained  (momentarily)  by  the  crystal 
can  be  deduced. 


A.  Finite  Strain  Fundamentals* 

Denote  the  coordinates  of  a  mass  element  in  an  initial  (unstrained) 
coordinate  system  by  a-j ,  and  the  coordinates  in  a  final  (strained)  system 
by  x.,-,  with  the  transformation  given  by, 

*i  =  x.j(t,  a^ ,  a£ ,  a3) ,  i  -  1,  2,  3 

where  (2.13) 

a;  =  W  ar  a2*  a3^ 


and  tg  is  a  reference  time.  The  x^  are  thus  Eulerian,  or  spatial,  coordinates 
and  the  a.  Lagrangian,  or  material,  coordinates. 


:h i s  transformation  one  can  derive  an  expression  for  the  ratio  of 


specific  volumes: 


(2.14) 


*This  section  is  a  summary  of  portions  of  the  theory  as  presented  by 
THURSTON  (46) . 
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J  is  thus  the  determinant  of  the  Jacobian  of  the  transformation,  or  the 
"functional  determinant." 

The  strain,  Njk,  is  defined,  somewhat  arbitrarily,  from  the  difference 
in  the  squares  of  the  lengths  of  line  elements  by: 


2N4b  da.  da,  =  dx.  dx.  -  da.  da. 
Jk  j  k  ii  ii 


(2.15) 


3X4  3X,* 

Njk=1'2 


Here  and  in  the  following  the  Einstein  summation  convention  for 
repeated  subscripts  applies.  6-^  is  the  Kronecker  delta. 

Expanding  the  internal  (strain)  energy  in  a  power  series  in  the 
strains,  or.e  obtains  (at  constant  entropy): 


P0[E(N,S)  -  E(0,S)]  =  1/2  C=jkl  N13  Nkl  +  1/6  c,Jkl„  N,  .  Nk,  N„ 


+  1/24  ‘ijktap,  Nij  Nkl  %  + 


(2,16) 


In  this  expression  the  c?^  ...»  represent  the  second  and  higher 
order  isentropic  elastic  stiffness  coefficients.  The  first-order  term  is 
missing  since  the  reference  state  is  considered  to  be  one  of  zero  stress  and 


strain. 


We  now  define  quantities,  called  thermodynamic  tensions,  by 


lij  = 


(2.17) 


In  terms  of  these  quantities  the  elastic  constants  are 
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and  similarly  for  the  higher  order  coefficients.  Consequently, 

p0  dE  =  tu  dNij  (dS  =  0) 

Finally,  the  equilibrium  (non-dissipative)  components  of  the  stress 
are  obtained  from  the  thermodynamic  tens:",  j  the  relations. 


3xk  axm 

aa7  S’, 

J  l 


(2.18) 


The  above  formulas  provide  isr- tropic  constitutive  relations  in  terms  of 
the  elastic  stiffness  coefficients.  Other  forms  u'  constitutive  relations 
can,  of  course,  be  derived  in  a  similar  fashion. 

Low  pressure  acoustic  measurements  yield  a  mixed  third-order  constant 


of  the  form: 


Cijkmpq 


...  /3cijknu 
“  '  3N  }T 

pq  t 


where  the  subscript  T  means  the  derivative  is  taken  at  constant  temperature. 

The  corresponding  purely  isentropic  constant  is  given  by: 

3  C5  *  1/ 

cijkmpq  “  Cijkmpq  +  ^T/,po  Ct^  ckmpq  auv^Cijkmrs  “rs  “  ^  iff  (2.19) 


where  Ct  is  the  specific  heat  at  constant  tension  and  the  auv  are  thermal 


expansion  coefficients, 


-  f8Wuvl 


3T  't 


In  view  of  the  symmetry  of  the  stress  and  strain  tensors,  the  number 
of  subscripts  can  be  reduced  by  adopting  the  following  convention: 


11  -  1 


32  ->•  4 


22  -v  2 


31  +  5 


33  3 


21+6 


This  convention  is  employed  in  the  following. 
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B.  Application  to  Uniaxial  Strain  in  Quartz 

We  assume  the  deformation  to  occur  in  the  X  direction  only.  The 
coordinate  transformation  is  accordingly, 

x1  =  (1  -  y)a- 


Fornvjl as  ( 2 . 14X..rOu*gn( 2.13)  then  gi ve : 


J  ■  V/'.Q  -  1  -  Y. 

N1  =  y(v/2  -  1) 

p0Ce  ~  ^  =  */2  cn  ^'i2  +  cm  +  ^24  cim  *  *  *  * 

tk  =  ci;,  N1  +  1/2  cllk  N-j2  +  1/6  clllk  N-j3  +  ...  (k  =  1,  2  ...  6) 
or  writing  out  the  components: 


^5  =  c15  ^■j  +  cii5  +  .  .  . 

t6  =  c16  N1  +  ^2  c116  +  *  *  * 
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The  stress  components  are  then: 

<*]  =  (1  -  v)t1 

<>2  =  O  -  y)-1  t2 
03  =  (1  -  y)'!  t3 


04  =  (1  -  y)-1 

o5  =  t5  (2-20) 

a6  =  l6 


For  alpha  cuartz  compressea  in  the  X-direction  the  above  formulas  are  correct 
as  they  stand.  For  compression  in  other  directions  the  proper  translation 
of  subscripts  must,  of  course,  be  made  to  indicate  the  correct  constants. 

The  above  formulas  have  been  applied  to  uniaxial  compression  of  X 
and  2-cut  cuartz,  using  the  second  and  third  order  constants  determined  by 
McSKIMIN,  £t  £l*(39)  and  THURSTON,  et  al_.  (40).  Values  of  these  con¬ 
stants  are  shc-wn  in  Table  HI. 

The  resulting  curves  are  plotted  inFigs.  2.6,  2.7  and  2.9.  The  values  of 
shock  velocity,  Us  and  particle  velocity,  up,  of  Fig.  2. 6  were  obtained  from 
the  Hugonict  relations: 


U 


s 


1/2 


and 


Up  -  [o(V0-V)]1/2 

The  predictions  are  seen  to  fall  outside  the  estimated  error  of  the 
shock  data,  indicating  that  the  fourth-order  term  contributes  significantly 

4»r\  A  AMA'AAW  /  AA  J  4*  Ia  /•  4*  ma  a  a  \  1  a  m  iiaI  iiaa  a  -C  xiiaA*  m  A"  tL*  ..«.!»  «  A  «.*  mUX 

w  cnc.yjr  ^anu  wic  scFcso;  au  laiyci  values  ui  suratn.  m  uiuuyti  I  u  iiiiyTn, 

be  thought  that  the  discrepancy  is  due  in  part  to  the  use  of  isentropic 
second-or.'-r  moduli  and  mixed  isentropic-isothermal  third-order  moduli  to 


predict  Hu^oniot  states,  for  which  internal  energy  is  greater  than  for 
isentropic  compression,  a  straightforward  calculation  shows  that  the  errors 


TABLE  III 

Elastic  Moduli  of  Quartz* 


Modulus 

Value 

(ion  dyne/cm^) 

Reference 

2nd-order 

cn** 

8.757 

(39) 

c12 

0.704 

U 

c13 

1.191 

II 

c14 

-1.804 

11 

c33 

10.575 

II 

3rd-order 

clll 

-21.0 

(40) 

c112 

-34.5 

II 

C1 1 3 

1.2 

II 

C1 1 4 

-16.3 

II 

C1 33 

-31.2 

II 

c333 

-81.5 

II 

4th-order 

cllll 

1705 

Present  Work 

c3333 

1849 

i; 

Tr.e  second-order  constants  are  isentropic,  the  third-order 
are  mixed  isothermal,  isentropic  constants,  and  the  fourth- 
order  are  Hugoniot  constants,  (see  text). 

The  c*| "|  constant  used  is  appropriate  for  open  circuit  com¬ 
pression,  i.e.,  at  constant  electric  displacement,  D. 
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tnus  produce-  are  negligible. 

Trie  c.fferences  between  the  purely  isentropic  third-order  moduli  and 
the  mixed  mc-uli  given  in  Table  inpan  be  calculated  from  tq.  (2.19) 

The  temperature  coefficients  of  expansion,  as  given  by  Il-vSON  (47) 

are: 

a3  =  7.8  x  10"6,  a-!  =  a2  =  14.3  x  10"6 

ar.d  the  expression,  due  to  Westrum,  reported  by  McSKIMIN  (39)  for  the 
specific  hea:  is: 

CP  (T)  =  Cp  {V  +  (T  '  VC1  +  {T  +  Tc)2c2  +  (T  '  Tc)3c3  a  *  *  * 

(77,4°K  <  T  <  298°K) 

where 

Tc  =  190°K 

Cp  (Tc)  =  5.189  x  106  erg/g°K 
C]  =  2.444  x  104  erg/g°K 
C£  =  -4.126  x  TCP  erg/g°K 
C3  =  5.327  X  10“2  erg/g°K 

taking 

T  =  25°C,  pQ  =  2.6485  g/cm3,  Cp  =  7.42  x  106  erg/g°K, 

■ 

and  estimating  (^—^-)  from  McSkimin's  data  taken  at  25°C  and  -195.8°C  to  be 

-  i 

of  the  order  of  -1  x  10®  dyne/cm3  °K  we  find  the  difference  given  by  Eq.  (2.19) 
for  the  c333  constant,  for  example,  to  be  of  the  order  of  5  x  1C3  dyne/cm3. 
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Tr.is  is  fc*r  orders  of  magnitude  less  than  C333.  Hence,  although  the  above 
calculation  is  hardly  accurate  for  the  cases  under  consideration,  only 
pathological  behavior  of  some  of  the  thermodynamic  variables,  a,  CQ  or 
{— £— )  ccjld  significantly  influence  the  results. 

The  difference  between  Hugoniot  and  isentropic  compressions  can  also 
be  shown  to  be  negligible.  For  compression  in  the  Z  direction  to  a  relative 
volume  of  3.9  the  strain  energy  given  by  Eq.(2.16)to  terms  of  third-order 
is  2.5  x  13^  erg/gm.  The  internal  energy  on  the  Hugoniot  is  2.8  x  109  erg/gm. 
Taking  Gruneisen's  ratio,  r,  to  be  approximately  1*,  the  stress  difference 
due  to  this  difference  in  thermal  energy  is  less  than  1  kbar--very  much  less 
than  the  c-served  stress  difference,  and  within  the  experimental  scatter. 

C.  Fourth-Order  Constants 

The  discrepancies  between  the  observed  data  and  the  predictions  based 
on  low  pressure  data  can  be  used  to  evaluate  fourth-order  coefficients.  This 
was  done  for  X  and  Z-cut  crystals  to  yield  the  values  of  c-j -j  1  -  and  c^-33  shown 
in  Table  13.  The  procedure  followed  was  to  fit  differences  between  the  data 
and  the  third-order  predictions  with  a  straight  line.  Because  of  the  large 
differences  in  pressure  range  and  experimental  precision,  this  method  proved 
to  give  an  adequate  fit  to  both  the  high  and  low  pressure  data.  No  adjust¬ 
ment  of  the  second  or  third-order  constants  was  necessary. 

The  fits  obtained  using  the  constants  up  to  fourth-order  are  shown  in 
Figs.  6,  7,  and  9.  Note  that  for  X-cut  crystals  the  slope  of  the  curve  in 
the  shock  velocity-particle- velocity  plane  (Fig/'S)  is  always  negative  when 
constants  only  up  to  and  including  third-order  are  used.  This  would  imply 

*AKDERS0N  (44)  gives  a  value  of  0.745  for  hydrostatic  compression.  Calcu¬ 
lated  for  tne  individual  components  with  the  assumption  Cp  =  Cv  gives 
I'll  —  T22  —  1.17;  T33  =  0.53. 
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a  shock  vjdVc  in  this  direction  is  unstable  end  spreads  as  ic  travels. 

With  the  a-dition  of  the  fourth-order  constant,  however,  the  slope  is  always 
slightly  positive.  Thus,  the  addition  of  the  fourth-order  term  results  ir.  a 
qualitative  difference  in  predicted  behavior. 

The  Us  -  Up  plot  for  Z-cut  crystals  is  nearly  a  straight  lira;  hcv.ever 
it  is  easi  j  shown  that  a  straight  line  does  not  accurately  fit  the  slope  at 
zero  particle  velocity.  Thus  the  straight  line  relation  often  assumed  in 
shock  studies  is  only  an  approximation  for  quart2  shocked  in  either  the  X  or 
Z  direction. 

It  is  also  easily  shown  that  the  Kurnaghan  form  of  equation  of  state, 
when  fittec  to  the  correct  slope  and  curv-ture  of  the  c  ~  V  curve,  (utilizing 
second  and  third-order  constants),  does  not  accurately  fit  the  higher  pressure 
data  and  is  therefore  an  approximation  only. 

These  statements  can  be  illustrated  by  examining  the  derivatives  of 
each  function.  Expanding  the  relation  for  o,  in  Eq.  2,20  in  terms  of  y 
yields: 

c  c 

a  =  Ci i  y[l  -  1/2  (3  +  -j.— — )y  +  1/6  (3  +  6  +  .  .  .]  (2.21) 

"  cn  cn 

A  linear  relation  between  shock  and  particle  velocity  of  the  form 

V  =  a  +  b  y, 
s  P 

ran  he  written,  by  means  of  the  Rankinp-Hugoniot  jump  conditions,  as 

a2  y  _ 

0  “  po  (I  -  by)2  » 

where  pQ  is  initial  density,  and  can  be  expanded  to  give 


a  -  pQa2  y[l  +  2by  +  3b2  y2  +  .  .  .] 


(2.22) 
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Finally,  the  one-dimensional-straif.  analogue  of  the  ‘vurnaghan  equation 

o  *  §  C(V0/V)B  -  1] 

can  be  expanced  to  give 

a  =  A  y[l  +  1/2  (3  +  1)Y  +  1/6  (B  +  1)(B  +  2)y2  t  .  .  .]  (2.23) 

Equating  the  derivatives  up  to  second-order,  we  have 

c-j -j  =  pQ  a2  =  A  (1st  order) 

-(3  +  -ail)  =  4b  =  (B  +  1)  (2nd  order) 

C11 

Evaluating  the  parameters  A,  b,  and  B  from  these  equations  gives 

X-cut:  A  =  pQa2  =  c^  =  8.68  x  10^ 

b  =  -0.15 
B  =  1.6 

Z-cut:  A  =  p  a2  =  c,,  =  10.575  x  lO^1 
o  o3 

b  =  1.177 
B  =  3.71 

With  these  values  all  three  expressions  have  the  same  slope  and  curvature  at 
zero  stress.  The  predicted  stresses  for  various  compressions  are  shown  in 
Table  IV. 

That  the  closed  form  expressions  are  approximate  is  hardly  surprising 
inasmuch  as  they  are  both  empirical  with  no  known  physical  basis.  Their 
value  is  that  they  both  are  two-parameter  functions  that  have  physically 
reasonable  shapes,  and  they  are  therefore  convenient  for  interpolation  and 
extrapolation  when  experimental  information  is  lacking. 
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2.34  Conclusions 

The  bnhavior  of  quartz  under  shock  loading  conditions  is  very  much 
different  from  that  of  metals,  as  was  pointed  out  by  Vlackerle.  The  elastic 
precursor  waves  are  an  order  of  magnitude  higher  and,  correspondingly,  so 
are  the  shear  stresses.  The  curve  labelled  X  in  Fig.  2.9isthe  stress  com¬ 
ponent  (base..  on  constants  to  third-order)  parallel  to  the  shock  front  when 
the  shock  propagates  in  the  2  direction.  The  maximum  stress  difference  is 
seen  to  exceed  100  kbar.  This  is  of  the  same  order  of  magnitude  as  the 
effective  shear  modulus;  consequently,  it  appears  that  quartz  momentarily 
exhibits  theoretical  yield  strength  under  dynamic  conditions. 

That  cohesion  of  the  material  is  destroyed  uoon  yielding  is  indicated 
by  the  close  agreemen  he  second  shocked  states  with  3ridgm-.n's  hydro¬ 
static  data.  There  is  no  indication  of  a  residual  shear  stress,  in  contrast 
to  the  case  for  metals  (49) . 

The  pronounced  stress  relaxation  shown  by  the  observed  variation  in 
amplitude  of  the  elastic  waves  and  the  apparent  dependence  on  tne  final 
pressure  is  quantitatively  larger  than  for  metals,  although  similar  quali¬ 
tatively. 

Evidently  shock  wave  methods  provide  a  valuable  supplement  to  low 
pressure  acoustic  measurements  in  determining  higher  order  elastic  constants, 
ac  least  for  ceramic  type  materials  which  sustain  large  amplitude  elastic 
waves.  Shoc.<  waves  are  inherently  more  suitable  for  higher  pressure 
measurements  than  are  acoustic  methods,  but  are  less  suitable  for  the  high 
precision,  low  pressure  measurements  required  to  evaluate  second-order  con¬ 
stants.  To  what  extent  shock  wave  techniques  are  capaole  of  measuring 
coefficients  other  than  the  principal  coefficients,  i.e.,  those  directions 
for  which  the  elastic  wave  is  purely  longitudinal,  requires  additional  study. 


III.  PHASE  TRANSITIONS 


3.1  Introduction 

Solid-solid  phase  transitions  are  classified  as  first 
order  if  a  volume  discontinuity  exists  at  constant  temperature 
and  as  second  or  higher  order  if  discontinuities  exist  only  in 
various  thermodynamic  derivatives.  We  are  concerned  here  only 
with  first  order  changes.  If  such  a  change  occurs  reversibly, 
the  transition  region  is  characterized  by  the  equality  of  the 
Gibbs  energies  of  the  two  phases.  This  leads  to  the  Clausius- 
Glapeyron  equation  which  relates  changes  in  pressure  and  tem¬ 
perature  in  the  mixed-phase  region: 

dp/dT  =  AS/Av  (3.1) 

The  consequences  of  this  relation  to  shock  wave  propa¬ 
gation  have  been  developed  in  some  detail  in  Reference  (7).  One 
of  the  important  results  is  a  relation  between  dp/dT  and  the 
slope  of  the  Hugoniot  in  the  coexistence  region  (8): 

(dp/dT)2  +  A(dp/dT)  +  B  -  0  (3.2) 

where 

A  =  2q-1/pm-31) 

B  Cpl/Tvl^M“0P 

0M  =  -  (l/v)(dv/dp)M  =  Hugoniot 

compressibility  in  mixed  phase 
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3]^  =  isothermal  compressibility 

cy^  =  volume  expansion  coefficient  at 
constant  pressure 

Cpi  =  specific  heat  at  constant  pressure 

Subscript  "l"  denotes  evaluation  in  phase  one 
at  the  phase  boundary. 

Measurements  made  on  bismuth,  iron  and  quartz  allow 
dp/dT  to  be  estimated  from  Eq.  (3.2).  Results  are  entered  in 
Table  V  where  static  values  are  also  given  for  comparison. 
Values  of  and  Cp^  at  atmospheric  pressure  and  room  tempera 
ture  were  used  in  the  estimates,  and  3-^  was  estimated  from 
shock  measurements. 


TABLE  V 

The  Coefficient  dp/dT  from  Static  and 
Shock  Wave  Experiments 


Type  of  Experiment 

Sample 

Bismuth 

Iron 

Quartz 

Static  measurements 

-.05  (9)* 

-.065  (LQ)* 

.018  01)* 

Shock  wave  measurements 
(calculated  from 

Eq.  (3.2)) 

-.067  *.045(7)* 

-.29  *.115  (7>v 

.225  (7)* 

*Pefer  to  sources  in  "Literature  Cited." 


From  the  table  it  is  apparent  that  shock  wave  data  for  iron  and 
quartz  do  not  agree  with  those  of  static  experiments.  The 
errors  in  shock  measurements  shown  in  the  table  are  estimated 
upper  bounds  based  on  rors  in  measurements  of  3^.  Therefore 
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it  is  very  difficult  to  attribute  the  discrepancies  for  iron  and 
quartz  to  experimental  errors  in  determining  the  slope  of  the 
Hugoniot  at  the  phase  boundary.  Conceivable  reasons 
for  the  discrepancies  are: 

1.  Inaccuracies  in  thermodynamic  data  (a^  and 

high  pressures.  Values  used  in  the  preceding  table 
are  those  at  atmospheric  pressures. 

2.  The  experimental  Hugoniot  curve  in  the  coexistence 
region  is  not  at  equilibrium.  There  may  be  some  non- 
equilibrium  rate-dependent  effects  influencing  the 
measurements. 

In  order  to  explain  the  discrepancies  of  Table  II, 
either  of  the  coefficients  a ^  or  Cp^  must  exhibit  a  two-  to 
four- fold  difference  between  its  value  at  high  pressure  and  that 
at  atmospheric  pressure.  Information  for  estimating  such  dif¬ 
ferences  is  not  generally  available.  Measurements  by  Bridgman 
below  20  kb  show  changes  of  not  more  than  20  per  cent  (12).  It 
seems  unlikely  that  changes  of  the  required  magnitude  will  occur 
at  these  higher  pressures. 

If  there  exist  rate  dependent  effects  in  phase  transition, 
we  should  be  able  to  observe  other  dynamic  evidence  in  wave  prop¬ 
agation.  For  example,  if  the  phase  transition  exhibits  a  relax¬ 
ation  process,  then  we  expect  a  decay  of  the  first  shock  wave 
with  travel  distance.  This  arises  because  the  delay  in  transi¬ 
tion  allows  a  mass  in  the  first  phase  to  momentarily  support  a 
higher  pressure  than  its  transition  pressure;  i.e.,  the  mass 
exists  in  the  extended  metastable  region  (Fig.  3.1).  This  is 
similar  to  stress  relaxation  of  an  elastic  precursor. 
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There  are  two  experimental  observations  of  the  decay  of 
the  first  wave  for  iron  (13,14). 

For  (v-quartz  Wackerle  (15)  did  not  observe  r  two-wave 
structure  corresponding  to  the  stishovite  transformation.  This 
may  indicate  that  the  transition  is  neither  rapid  nor  complete. 

In  order  to  study  transient  effects  in  shock  wave  prop¬ 
agation,  we  must  solve  the  flow  equations.  In  Lagrange  coordi¬ 
nates  for  plane,  one -dimensional,  time- dependent  flow  these  are: 

1.  Newton's  second  law 

pQ(Wat)h  +  (d(p+q)/dh)t  =  3  (3.3) 

where  h  =  Lagrangian  space  coordinate 
t  =  time 

p0  =  initial  density  =  l/v0 
q  =  viscous  stress 

Introduction  of  q  makes  the  flow  continuous  through 
the  shock  front  and  the  jump  conditions  unnecessary. 

2.  Continuity  equation 

p0(av/at)  -  au/ah  =  0  (3.4) 

3.  Energy  conservation 

eE/at  +  (p+q)  (av/at)  =■  0  (3.5) 

4.  Constitutive  relation 

In  order  to  solve  the  above  differential  equations  in 
tour  independent  variables  (p,  v,  u,  e) >  constitutive 
relations  are  required.  Ignoring  solid  rigidity,  which 
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should  play  a  secondary  role  in  the  problems  considered 
here,  these  will  consist  of  an  equation  of  state,  p  = 
P(v,e)>  an  expression  for  q,  and,  in  the  mixed  phase, 
a  statement  of  the  kinetics  of  transition. 

3.2  General  Constitutive  Relations 
for  Two- Phase  Flow 

We  assume  the  following  conditions  to  apply  in  a  small 
mass  element  when  both  phases  coexist: 

1.  Equal  pressure  exists  in  both  phases. 

2.  Temperatures  of  the  two  phases  are  equal ,  i.e.,  the 
heat  released  by  the  transition  is  instantly  re¬ 
distributed. 

3.  Particle  velocities  are  the  same  for  each  phase. 

4.  No  surface  energy  is  associated  with  the  interface 
between  the  two  phases. 

5.  No  heat  is  transferred  between  mass  elements  by 
conduction. 

Conditions  1  and  3  insure  that  Eqs.  (3.3)-(3.5)  apply 
for  the  coexistence  region  as  well  as  for  a  single  phase.  The 
only  change  needed  is  a  reinterpretation  of  v  and  E .  This  can 
be  done  as  follows: 

Regardless  of  phase  transition,  the  total  mass  of  a 
given  Lagrangian  volume  is  conserved.  Therefore,  if 
we  denote  by  v  the  total  specific  volume  of  the  mix¬ 
ture  of  two  phases,  then  from  condition  3  above,  v 
satisfies  Eq.  (3.4)  on  p.  3.4.  But  the  variable  v  must 
also  satisfy  the  following  relation  in  the  coexistence 

=  (1-or)  v^(p ,T)  +  o'V2(p,T) 


region: 

v(p,T) 


(3.6) 
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where  a  is  defined  as  the  mass  fraction  of  the  second 
phase  and  v^  is  specific  volume  of  the  i-th  phase  at 
p  and  T. 

Conditions  4  and  5  imply  that  the  total  internal  energy 
of  a  mass  element  is  the  sum  of  the  internal  energies  of  the  two 
phases : 

E  -  (l-o)E1  +  aE2  (3.7) 

where  is  specific  internal  energy  of  the  i-th  phase 
at  p  and  T. 

In  order  to  obtain  a  suitable  constitutive  relation  when 
the  phase  transition  has  a  finite  reaction  rate,  we  assume  first 
that  the  p,v,E  surfaces  in  each  phase  can  be  extended  smoothly 
into  metastable  regions  overlapping  the  equilibrium  region  of 
mixed  phase,  as  in  Fig.  3.1. 

Then  we  have  the  following  functional  forms: 

V1  =  V1<P>T) 

v2  =  v2(p,T) 

E-,  =  EAp,?) 

i.  y  i 

E2  =  E2(p,T).  (3.8) 

Then  from  Eqs.  (3.6)  and  (3.7)  we  have 

dv  =  (l-a)dv^  +  adv2  +  (v2“V^)dff  (3.9) 

dE  -  (l-ff)dE].  +  ordE2  +  (E2“E^)da’.  (3.10) 

From  Eq.  (3.8) : 

dvj^  =  (B»vi/?>T)pdT  +  (dv^Ap^dp  (3.11) 

i  =  1  or  2. 
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With  E  a  function  of  T  and  p  we  have: 
dEj.  =  (^Ei/^DpdT  +  (^EiAp)Tdp 

=  (C  £  -  p(avtAT)  )dT  -  (T(»Vi/aT)  +  p(avi/>p)T)dp 

(3.12) 

i  =  1  or  2 

where  =  specific  heat  of  i-th  phase  at  constant  pressure. 
Substituting  (3.11)  and  (3.12)  into  (3.9)  and  (3.10),  we  get 

dv  =  i^dp  4-  m^dT  +  n^d(y  (3.13) 

dE  =  j^dp  +  m2^  +  n2c^Q'  (3.14) 

where 


1 1  13  (1-or)  (Sv^/^p)<ji  +  Qr(ftV2/’ (3.15) 

m1  =  (1-Q-)  (^v3/^T)p  +  ff(av2/aT)p  (3.16) 

n^  =  v2  "  vi  (3.17) 

12  =  -(l-o-)  (T^Vj/aT^  +  p(avj/ap)T) 

-  a  (T(av2/*T)p  +  p(*v2/*p)T)  (3.18) 

m2  -  (1-or)  (C  ,  -  p(^v1AT)  )  +  o(C  2  -  p(*v2/aT)  ) 

(3.19) 

n2  -  E2  -  Ei  (3.20) 

and 

dS  =  -(p-kj)dv  from  Eq.  (3.5).  (3.21) 

Therefore,  in  principle,  (3.13)  and  (3.14)  can  be  solved  for  dp 
and  dT  if  d«y  is  given. 

Now  we  assume  the  following  relaxation 
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relation  whose  derivation  is  the  subject  of  the  next  section, 

do-  =  f(v,T,a)dt.  (3.22) 

Then  substituting  Eqs.  (3.21)  and  (3.22)  into  (3.13)  and  (3.14) 
and  solving  for  dp  and  dT,  we  get 
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and 

Cpi  =  specific  heat  capacity  at  zero  pressure. 

In  summary,  the  generalized  constitutive  relations  now 
consist  of: 

1.  Equations  of  state 

vi  5=5  vi(P»T) 

2.  Specific  heat  capacities 

cPt  ■  VP>T) 

3.  Zero-point  energy  difference 

„oo  oo 
Ei  -  E2 

4.  Relaxation  relation 

dcr/dt  =  f(v,T,ty) 

5.  Artificial  viscosity,  q. 

* ,  '  4 

To  complete  the  above  description  we  must  find  a  relaxation 
relation  and  an  expression  for  q. 


3.3  Irreversible  Thermodynamics  and 
Phase  Relaxation 

This  section  is  concerned  with  the  mechanism  of  phase 
transformations  in  the  solid  state.  It  should  be  possible  to 
describe  the  mechanism  of  phase  change  in  solids  in  terms  o£- 
interatomic  forces  by  use  of  kinetic  theory.  Actually,  because 
of  the  problem's  complexity,  no  such  quantitative  description 
has  been  achieved.  However,  very  successful  phenomenological 
theories  for  the  kinetics  of  phase,  change  have  been  devel¬ 
oped  (16,17,13,19).  Such  theories  are  of  two  kinds,  known  as 
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nucleation  and  growth  processes  and  martensitic  trans tor mat ions 
This  nomenclature  is  rather  unfortunate,  as  pointed  out  by 
Christian, because  growth  from  nuclei  also  occurs  during  mar¬ 
tensitic  reactions. 

In  nucleation  and  growth  processes  a  new  phase  grows 
from  critical  nuclei  at  the  expense  of  the  old  phase.  The  re¬ 
action  proceeds  to  completion  by  a  slow  migration  of  the  inter 
phase  boundary,  the  velocity  of  which  varies  markedly  with  tem¬ 
perature.  Most  atoms  have  different  neighbors  in  a  new  phase. 
On  the  other  hand  in  martensitic  transformation  the  reaction 
takes  place  by  coordinated  atom  movements  (e.g.,  shear-like), 
and  atoms  have  the  same  neighbors  after  transformation.  The 
latter  is  often  observed  in  rapid  cooling  of  alloys. 

The  basic  model  for  describing  nucleation  and  growth 
phenomena  is  based  on  the  following  conditions: 

1.  Steady  state  reaction 

2.  Isothermal  and  isobaric  process 

3.  Thermal  fluctuations  as  driving  force 

4.  Boltzraan's  law  for  the  relative  probability  of 
different  energy  states. 

The  expression  of  this  theory  is: 

1.  Nucleation 

I  (nucleation  frequency) a  exp(-AG*/kT) 
where  AG*  is  defined  to  be  the  critical  net  free 
energy  change  on  forming  the  nucleus  of  the  new 
phase. 
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2 .  Growth 

There  are  many  expressions  for  the  growth  law.  An 

often  used  form  is: 

x  =  1  -  exp(-ktn) 

where  x  =  the  volume  fraction  transformed 

k  =  a  complicated  function  of  the  nucleation 
frequency  I 

t  =  time 

n  =  constant  of  magnitude  3  to  4. 

The  principal  difficulties  in  applying  this  model  to  the  shock- 
induced  transition  are  the  conditions  imposed  in  deriving  the 
rate  equation.  In  a  shock  front  the  situation  is  neither  iso¬ 
thermal  nor  isobaric,  nor  is  it  steady  state. 

A  description  of  the  nucleation  of  martensitic  trans¬ 
formation  is  given  (20)  for  an  athermal  case  (e.g.,  quenching), 
but  the  result  is  less  conclusive  than  that  for  a  nucleation 
and  growth  process.  This  transformation  has  a  negligible  free 
energy  barrier  (21)  and  proceeds  very  rapidly  to  completion. 
Therefore  at  present  the  quantitative  information  obtained  is 
mainly  concerned  with  the  amount  of  transformation  with  respect 
to  the  cooling  rate  but  not  with  time. 

It  is  tempting  to  use  the  martensitic  transformation  in 
describing  shock-induced  transformations  because  of  its  high 
speed.  There  are  some  who  suggest  the  martensitic  transition 
for  iron  (a  -*  e)  in  shocks  (22.23).  However,  a  quantitative 
description  of  the  rate  equation  is  not  yet  available  for  inclu¬ 
sion  in  the  constitutive  relations. 
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As  more  and  better  data  on  shock- induced  transition  be¬ 
come  available,  it  will  become  important  to  develop  some  of  the 
above  models  to  describe  the  kinetics  of  phase  transformation  in 
shock  waves.  At  present  it  is  barely  established  that  rate 
effects  in  shock-induced  transformation  exist,  and  we  seek  no 
.ore  than  a  qualitative  description  of  the  effects  of  kinetics 
on  the  wave  structure  and  some  rough  numbers  for  the  magnitude 
of  the  reaction  rate  involved  in  the  observed  processes. 

To  this  end  we  select  a  quite  different  approach,  which 
is  in  some  sense  better  founded,  though  formally  limited  to 
small  deviations  from  equilibrium,  in  irreversible  thermodyna¬ 
mics.  In  this  approach  there  exist  none  of  the  conditions  im¬ 
posed  in  the  above  models,  and  it  suggests  a  simple  relaxation 
law  for  the  transformation. 

dcy/dt  =  (or  ^  -  a)/r  (3.26) 

where 

eq 

a  =  equilibrium  value  of  <y,  or  a  completely 
relaxed  value 

t  15  relaxation  time. 

The  derivation  of  the  above  relation  from  irreversible  thermo¬ 
dynamics  follows. 

When  a  solid- solid  transition  occurs  reversibly  (or  in 
equilibrium)  by  adiabatic  compression,  the  entropy  of  a  given 
system  is  constant  and  the  process  is  called  isentropic.  If  the 
transition  is  not  reversible,  the  entropy  will  increase.  How¬ 
ever,  ordinary  thermodynamics  does  not  give  the  precise  amount 
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of  increase;  rather,  it  tells  only  the  direction  of  increase: 

ds  sO 

irreversible 

Basically,  the  idea  of  irreversible  thermodynamics  is  to 
replace  the  inequality  by  an  equality  so  that  we  can  determine 
the  increase  of  entropy  caused  by  irreversible  processes. 

To  establish  the  desired  equality  we  assume  that,  al¬ 
though  the  total  system  is  not  in  equilibrium,  there  exists 
within  a  small  mass  element  a  state  of  local  equilibrium  for 
which  the  total  entropy  change  per  unit  mass,  ds,  is  expressed 
by  the  Gibb's  relation  (24). 


=  mass  fraction  of  the  i-th  component. 

This  means,  for  the  case  of  a  single-component  phase  transition, 
that  when  two  phases  are  not  in  equilibrium  with  one  another,  we 
interpret  them  as  if  they  are  two  different  components  with  dif¬ 
ferent  specific  Gibb 5 s  functions  and  each  is  in  local  equilibrium 
satisfying  Eq.  (3.27). 

There  is  an  additional  constraint  on  Eq.  (3.27)  for  a 
single-component  system  from  total  mass  conservation.  Since 
the  total  mass  is  constant: 


da^  +  da?2  ~  0* 


(3.28) 


J 
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We  redefine  a 2  =  =  1  "  a  and  use  o-  hereafter  as  the  sole 

reaction  variable.  Then  from  Eqs.  (3.27)  and  (3.28)  we  get: 

Tds  =  de  +  pdv  -  (3.29) 

So  far  we  have  been  looking  at  an  isolated  element  of 
mass.  But  if  we  assume  that  this  mass  element  is  part  of  a 
continuous  medium  which  is  in  a  state  of  non-uniform  stress, 
strain  and  motion,  then  the  variables  e  and  v  must  satisfy  the 
energy  and  mass  conservation  equations,  (3.4)  and  (3.5).  Sub¬ 
stituting  these  equations  into  Eq.  (3.29)  we  get: 

T(ds/dt)irrev>  =  -vq(au/*x)  “  (p.2“M>x) dof/ dt  (3.30) 

We  assume  no  heat  to  be  deposited  from  the  outside,  so  the 
entropy  change  is  entirely  due  to  the  internal  irreversible 
process.  It  should  be  noted  that  the  Lagrangian  derivative  used 
in  Eqs.  (3.4)  and  (3.5)  is  identical  to  the  convective  deriva¬ 
tive  implied  in  Eq.  (3.30).  If  we  look  closely  at  Eq.  (3.30), 
we  can  see  the  sources  of  irreversibility.  These  are  a  chemical 
affinity  (y. 2“W- 1)  >  (25),  and  a  velocity  gradient.  These  quan¬ 
tities  are  called  "forces"  in  irreversible  thermodynamics  and  are 
denoted  by  X^(i=l,2. . .) .  The  phenomena  caused  by  these  forces, 
such  as  phase  changes,  are  called  "fluxes"  and  are  described  by 
J^(i-1,2. . .) .  Then  Eq.  (3.30)  becomes: 

2 

T<ds/dT)irrev.  =  l  J ft  (3.31) 

i=l 


xx  =  -au/*x 

x2  =  *'GJ-2~fAl) 


where 


(3.32) 

(3.33) 
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Furthermore,  from  Eqs.  (3.36)  and  (3.37)  we  get: 

T(ds/dt)irreVa  =  g11X1  +  (gi2+g21^XlX2  +  g22X2‘ 

The  second  law  of  thermodynamics  requires  that  the  left-hand 
side  is  always  positive  or  zero,  so  we  must  have  the  following 
relations  for  the  constants  g^j : 

gll  s  °>  g22  s  0 

2  ^Sng22  *  lg12+g2l!* 


But  this  is  all  we  get  from  thermodynamics;  it  will  not  give  us 
the  magnitudes  of  the  coefficients,  because  they  are  character¬ 
istic  numbers  which  depend  upon  materials.  A  complete  under¬ 
standing  of  the  reaction  mechanisms  would  provide  a  foundation 
for  calculating  the  g  „ .  For  the  present  we  assume  the  reaction 
rate  to  be  independent  of  shear  stress  and  set  g-^  =  0*  &n  and 
g22  are  chosen  from  available  data  and  for  computational  con¬ 
venience,  respectively.  Then  Eqs.  (3.38)  and  (3.39)  reduce  to: 


dor/ dt  = 


(3.40) 


q  =  g22(i/v)(3u/ax) .  (3.41) 

Since  we  know  the  chemical  affinity  (pt, 2“M- 1)  is  zero  at 

eq 

the  equilibrium  state,  we  can  find  the  equilibrium  value  a  of 
a  as  a  function  of  two  independent  thermodynamic  parameters. 

For  example: 

*eq  - 


T). 


Then  for  a  small  deviation  from  equilibrium,  we  can  expand  the 
affinity  in  a  Taylor's  series: 


^2“^!  =  (u2"Hx)eq  + 


a(n  -p,  )eq 


v  ,T 


(ry-o»eq) 


Neglecting  the  higher-order  terms: 


>(n2-u1)®q 

P-  9~U  -i  ~  ( - ) 

1  L  v,T 


(o"-o'eq)  . 


(3.42) 


Equations  (3.40)  and  (3.42)  give: 

dry/dt  =  (<y-cyeq)/T  (3.42) 


Uhere  t  is  a  new  constant  to  be  determined  from  comparisons  of 
calculations  with  experiments. 

Eq.  (3.41)  shows  that  q  has  the  character  of  a  viscous 
force,  as  stated  earlier.  It  is  necessarily  proportional  to 
the  first  power  of  the  strain  rate  because  of  the  use  of  the 
linear  phenomenological  law.  This  first-power  dependence  is 
retained  in  the  computations,  but  since  its  primary  purpose  is 
to  smooth  the  shock  transition  it  will  be  artificially  modified 
after  the  manner  described  by  Richtmyer  and  von  Neumann  (27): 
The  coefficient  g  ^  made  proportional  to  the  cell  thickness 
used  in  the  computation. 

Then 


where  CL 


q  =  C^(Ax/v)  j*u/*xj 


(3.44) 


«  constant. 
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This  completes  our  general  description  of  non- equilibrium 
two-phase  flow.  In  the  following  section  we  approximate  the 
constitutive  relations  in  order  to  apply  them  in  numerical  cal¬ 
culations  . 


( 


IV.  CONSTITUTIVE  RELATIONS  FOR  IRON 

4.1  Approximations  for  Two-Phase  Flow 

Three  approximations  were  made  in  order  to  simplify  the 
application  of  the  principles  of  Section  III  to  the  m  -  z  tran¬ 
sition  in  iron.  First,  the  equilibrium  volume  change  at  con¬ 
stant  p  was  assumed  constant  in  the  coexistence  region.  This 
approximation  is  supported  by  compressibility  measurements  at 
room  temperature  which  result  in  the  following  values  (10,28): 

&2  (at  192  Kb)  =  4.94  x  lO^O-Jb)”1 

(at  132  Kb)  =  5.1  x  10"1(Mb)'1 

where  9^  is  the  compressibility  of  the  i-th  phase. 

Then 


V2(P >T)  -  v^(p,T)  =  -.004  cc/g  at  every  p  and  T.  (4.1) 

Then  from  Eqs.  (3.6),  (3.26)  and  (4.1)  we  get: 

dv^  =  dv  -  (v2”V^)do'  (4.2) 

d<y  =  (ote<^-a)dt/r  (4.3) 

Then  v-^  can  be  calculated  from  the  above  equations  provided 
v  and  c*e<:*  are  known.  In  the  computation  process  v  is  given  by 
the  continuity  equation  and  <*ec^  by  Eq.  (3.6). 
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aQq  =  ((v-v1)/(v2-v1))eq  =  (v-v®q(T))/(v2-v1)eq  (4.4) 

Since  (v2~v^)e<^  is  constant,  the  temperature  dependence  of  <*eq 
is  in  v^(T)  .  Our  second  approximation  is  to  replace  v-^(T)  by 
the  value  v^(TQ)  at  room  temperature  (T  ) .  Since  the  tempera¬ 
ture  is  only  slightly  greater  than  room  temperature,  the  error 
involved  in  this  approximation  is  negligible  except  in  the 
immediate  vicinity  of  the  transition  pressure,  130  Kb. 

Since  v^  can  be  determined  from  Eqs.  (4.2)  and  (4.3)  we 
can  use  the  form  p  =  p(v^,T)  to  calculate  pressure  if  we  know 
the  temperature.  The  third  approximation  we  made  is  the  tem¬ 
perature  calculation. 

If  the  phase  transition  is  made  at  equilibrium,  pressure 
and  temperature  in  the  coexistence  region  must  satisfy  the  rela¬ 
tion: 


dp/dT  =  constant  at  fixed  p  or  T. 

Over  a  wide  range  of  temperatures  the  coefficient  is  practically 
constant  (10).  But  if  we  assume  that  the  coefficient  dp/dT  is 
constant  even  in  the  non-equilibrium  coexistence  region,  we 
shall  see  presently  that  the  energy  conservation  law  alone  is 
then  sufficient  to  determine  the  temperature  change.  Since  E  is 
given  by  Eq.  (3.7)  the  internal  energy  is  a  function  of  p,  X  and 
a.  On  replacing  a  by  Eq.  (3.6),  it  is  a  function  of  p,  T  and  v. 
But  if  we  assume  that  dp/dT  is  constant,  p  and  T  are  no  longer 
independent  and  E  is  a  function  of  T  and  v  only.  The 
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error  resulting  from  this  approximation  has  not  yet  been 
evaluated.  It  should  be  small  since  AT  is  small.  From 
this  assumption  we  relate  changes  in  temperature  to  changes  in 
specific  volume  as  follows: 

dE  »  (6E/&T)vdT  +  (3E/*v)Tdv 

Substituting  into  this  relation  from  Eq.  (3.21),  we  find  in  the 
mixed  phase  region: 

(6E/dT)vdT  “  -(p+q)dv  -  (SE/6v)Tdv. 

Defining  (SE/3T)  as  C  and  using  the-  following  identity, 

V  V  jID 

T(aP/&T)v  -  p  +  (3E/3T)v  , 

this  becomes 

dT  *  -  (T(3p/*T)V  +  q)dv/Cv>m  ,  (4.4) 

The  calculation  of  Cv  m  in  the  coexistence  region  is  given  in 
Appendix  I.  It  should  be  noted  that  in  the  coexistence  region 
(Sp/ST)v  is  equal  to  dp/dT  because  of  the  assumed  dependence  of 
p  and  T. 

Our  assumption  or  internal  energy  dependence  (E(v,T)) 
is  certainly  true  for  a  single  phase,  sc  the  form  of  Eq.  (4.4) 
is  valid,  and  if  wc  know  the  values  of  C  and  (5p/dT)__  in 
either  a  single  phase  or  a  mixed  phase,  the  temperature  cal¬ 
culation  is  only  a  matter  of  substituting  different  values  of 
these  parameters,  depending  on  the  phase  region. 
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4.2  Review  of  Experimental  Information 

There  are  three  main  static  experiments  on  the  equation 
of  state  (12,10,29)  and  the  isothermal  p-v  relations  are  in 
close  agreement,  within  the  range  ot  experimental  errors,  with 
shock  wave  measurements  ot  the  Hugoniot  at  low  pressures.  Ve 
use  the  equations  of  state  determined  from  shock  experiments. 

There  are  tour  major  measurements  on  shock  wave  propaga¬ 
tion  in  iron.  The  one  by  Bancroft  et  al.  (28),  is  on  the 
Hugoniot  pressure-volume  relations  above  the  transition  pressure 
and  it  reveals  a  discrepancy  in  the  magnitude  ot  the  gradient 
dp/dT  compared  with  the  static  value  (/).  The  other  three  are 
concerned  with  transient  effects  and  the  thickness  of  the 
second  shock: 

1.  Bancroft  et  al.  observed  the  two-wave  structure  and 
determined  the  pressure-volume  relation  in  the  second 
shock  by  jump  conditions. 

2.  Smith  measured  the  thickness  (.02  mm)  of  the  second 
shock  front  in  recovery  experiments  by  hardness  methods 
and  concluded  that  the  duration  of  the  transition  zone 
is  in  the  order  of  .001  ^sec  (30), 

3.  Novikov  et  al.  observed  the  two-wave  structure  in  iron 
by  a  capacitance  method  (14) .  They  concluded,  from  the 
rise  time  of  the  second  shock,  that  the  duration  of  the 
transition  region  is  .2  -  .3  psec.  This  duration  is 
not  exactly  equal  to  the  relaxation  time  t,  but  as  it 
is  essentially  governed  by  the  relaxation  of  the  phase 
transition,  it  should  be  of  the  same  order  of  magni¬ 
tude.  They  also  observed  a  pressure  drop  behind  the 
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first  shock  for  thinner  samples  (no  information  is 
given  about  the  thickness) ,  and  explain  that  it  is  due 
to  relaxation  processes. 

4.  Minshall  observed,  by  a  pin  technique,  decay  of  the 

first  shock  over  a  distance  of  5  cm  (.8  cm  *  5.8  cm)  (13). 
No  estimation  of  transformation  time  is  given,  nor  is 
the  explanation  of  the  decay.  However,  if  we  assume 
this  decay  to  be  due  to  the  relaxation,  the  order  of 
relaxation  time  must  be  about  5-10  p,sec. 

All  of  the  above  reports  agree  about  the  existence  of 
the  transition,  but  as  far  as  the  relaxation  time  is  concerned, 
they  make  no  suggestion  of  a  particular  value  to  be  used  in  the 
calculation. 

In  numerical  procedures  we  can  use  any  relaxation  time 
to  study  the  effect  of  phase  change  on  shock  wave  propagation, 
but  we  made  an  arbitrary  choice  of  1/3  psec  for  most  of  the 
calculations,  based  on  consideration  of  the  experiments  by 
Novikov.  In  the  study  of  the  decaying  precursor  we  used  three 
relaxation  times,  .1,  1/3,  and  1  psec. 


4.3  Equation  of  State  of  Iron 

The  equation  of  state  of  the  first  phase  is  taken  to  be: 


p(v^,T)  ~  a^(M^-l)+a2(T|^-l)^r-32(T|^‘-l)^'-Cv^(T-T0)r/v^ 


/  A  C\ 

v-t.-v 


where 

C  1  =  specific  heat  at  constant  volume,  phase  1, 
vi  assumed  constant 

T  =  some  temperature  above  which  CL-,  is  constant, 
taken  as  room  temperature  here. 
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a^  =  constants 

T  =  Gruneisen  function,  assumed  constant 

Hi  =  v0/v1. 

This  is  a  form  of  the  Mie-Gruneisen  equation  used  by  Al'tshuler 
g.tL.a.l»  (31). 

The  coefficients  a^  can  be  determined  from  the  poly¬ 
nomial  fits  of  the  Hugoniot  curve  (1)  or  from  static  measure¬ 
ments.  However,  for  the  case  of  iron  there  is  no  appreciable 
difference,  below  200  Kb,  between  the  isotherm  and  the  Hugoniot 
centered  at  room  temperature.  For  example,  the  temperature  rise 
along  the  Hugoniot  from  0  to  130  Kb  is  20°C  (23)  ,  which  contrib¬ 
utes  only  about  1.3  Kb  to  the  total  pressure.  This  difference 
is  less  than  experimental  error  for  static  and  shock  measure¬ 
ments  in  general  (29).  Therefore,  we  can  substitute  the  Hugoniot 
as  a  room  temperature  isotherm  in  the  equation  of  state.  These 
and  other  equation  of  state  parameters  are  given  in  Table  VI. 

The  values  of  a^  listed  in  Table  VI  are  determined  from  the 
least  square  fit  of  existing  data  below  130  Kb  (32).  Since 
errors  in  the  experiments  are  larger  than  the  thermal  pressure, 
this  will  not  give  any  inconsistency  in  the  equation  of  state. 
When  we  speak  of  the  temperature- independent  equation  of  state, 
we  mean  the  isotherm  at  T0. 

Fig.  4.1  shows  the  isotherm  at  T0  in  terms  of  relative 
volume.  Bridgman's  data  at  room  temperature,  extrapolated  to 
high  pressures,  are  drawn  for  comparison.  The  difference  at 
high  pressures  is  mainly  due  to  inaccuracies  encountered  in 
extending  Bridgman's  data  to  such  high  pressures.  Since  the 


volume  change  is  constant  at  fixed  p  and  T,  the  isotherm  of 
the  second  phase  (e  -  phase)  can  be  easily  found  by  shifting 
the  first  phase  by  the  amount  of  the  volume  change  (v2_v^) . 


TABLE  VI 


Physical  Data  for  o>  -  iron 


Parameter 

Values 

Dimension 

Reference 

vG  (initial  volume) 

.1275 

cc/g 

(33) 

a-,  (thermal  expansion 
coefficient) 

36.3  x  10"6 

i 

l/K°(degK)"1 

(33) 

Cvi  (heat  capacity) 

.4447x10" 5 

Mbcc/g°K 

(33) 

p^  (transition  pressure) 

.130 

Mb 

(10) 

(dp/dT)ra  (equilibrium) 

-.000065 

Mb/K° 

(10) 

Av  (volume  difference) 

-.004 

cc/g 

(10) 

al 

1.667 

Mb 

(32) 

a2 

3.4 

Mb 

(32) 

a3 

0 

Mb 

(32) 

r 

1.6 

•  •  • 

(34) 

Cv  ,m 

.46xl0"5 

Mbcc/g°K 

* 

To 

300°K 

°K 

*See  Appendix  III. 

Once  the  equation  of  state  is  known,  the  expression  for 
eq 

o'  can  be  easily  found  from  Eq.  (4.4),  which  includes  the  room 
temperature  approximation  for  v^(T).  Suppose  we  specify  the 
room  temperature  transformation  by  the  isotherm  AB  in  Fig.  4.1; 


Relative  Volume,  v/v 
’  '  o 


Fig.  4. 1.  — Temperature  Independent 
Equation  of  State  of  Iron  and  o<ecl 


I 


69 

GCJ 

a  is  given  by  the  relations: 


= 

0 

if 

V  >  V 

A 

aeq  = 

(v-vA)/<vB-vA) 

if 

VB  <  V  <  VA 

(4.6) 

1 

if 

v  <  V 

The  graph  of  (YecJ 

is  given  in  Fig. 

4.1. 

As  seen  in  Table  VI,  we  assume  the  constancy  of  physical 
data,  such  as  0  ,  r  and  so  on,  regardless  of  pressure.  We  use 
the  equilibrium  value  (-.065  Kb/°K)  for  dp/dT  in  the  coexistence 
region  unless  otherwise  stated. 


V.  SHOCK  PROPAGATION  IN  IRON 


j.l  Difference  Equations 

The  system  of  difference  equations  used  here  is  based 
on  one  described  by  Wilkins  (35)  ,  with  the  yield  stress  set  to 
zero.  Space  is  divided  into  points  and  cells  as  in  Fig.  5.1. 

The  particle  velocity  and  the  current  position  of  the  Lagrangian 
coordinate  are  defined  on  the  points  with  integer  label  ) ,  and 
other  variables  are  assigned  to  the  cells  with  half-integer 
labels  j+%,  etc. 


Point  Point 


j+i 

Cell 

Cell 

x(j)  x(j+l) 

Fig.  5 .1.- -Difference  Scheme  in  Space 

Time  differences  are  also  staggered.  Particle  velocity 
and  q  are  defined  at  half- integer  times,  n4%,  and  other  vari¬ 
ables  are  defined  at  integer  times,  n. 

The  computational  sequence  for  general  interior  points 
and  cells  is  given  in  Fig.  5.2. 
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0  New  calculations 


•  Old  calculations 


First  step 


\  /  v 

\  /  \ 


/ks 

/  f  \\  v  \  \ 

//  \  \p//  \\ 

• - - - 1 . W _ L _ \ _ 


V 

q 


Second  step 
(volume  only) 

Third  step 


C  Fourth  step  (p,T) 


J-f  j  J+f  3+1 

Fig.5.2. — Computational  Sequence 


The  flow  equations  and  the  constitutive  relations  are  approx¬ 
imated  according  to  the  above  difference  scheme.  The  difference 
equations  are  given  in  the  order  of  computation.  The  subscript 
or  superscript  "o"  refers  to  the  initial  state  (p=0)  of  the 
first  phase. 

1.  Equation  of  motion,  Eq.  (3.3): 


n+%  n-%  ,  At  n  v. 

u3  U1  '  s?  (  j+*>  ' 

J 


where 


n  ,  n  ,  n-%. 

"(P  +q 


(5.1) 


(5.2) 


=  total  stress  in  the  x- direction. 


n  „n 


n  n 


xV+1-xV  x.-x.  x 

wHr^  +  cV-)i 


=  average  mass  at  j 


At  the  left  boundary 


(5.3) 


(5.4) 


The  new  coordinate  is  given  by: 

xn+l  =  xn  +  Un-Hs  4t. 

J  J  J 


(5.5) 


2.  Continuity  equation: 


v?S  -  v^+4tO^<-uf)  (5-6) 


where 


=  p? , v(x° , n -x?)  =  mass  in  the  cell  j+%,  (5.7) 

H3*^  3+1  3 


po  = 


=  initial  density. 


Linear  viscosity: 


n-t%  o  n+%  .  n+% 

<Jj4%  “  CL  Pj+%  ^3+%  1Uj+l 


un-^  <  un-^ 
n+%,  - _ f  3+1  3 


for( 


(5.8) 


0  otherwise. 


Here 


-  2vo/(vj+% + vW 


(5.9) 
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4.  Constitutive  relations: 

The  relaxation  equation: 

a^tv  =  a\?,i  +  (” - ~)  ,At.  (5.10) 

3-r%  J+3  v  T 

t  is  the  characteristic  relaxation  time  and 
is  assumed  to  be  constant. 

The  specific  volume  of  the  first  phase  is: 


„n+l  _  ..n+1  /  „  \  n+1 

vi,j+*  “  V.i-Hs  ' 


(5.11) 


Temperature  calculation: 


‘  T"^  +  CCT]^(v^-v;+%) 

- 

where  C  =  -(l/C  )(?>p/AT)v  (convenient  for  a 

mixed  phase) 

or  C  =  -r/v  (convenient  for  a  single  phase) 
and  r  is  the  Gruneisen  coefficient.  The  value 


(5.12) 


of  Cv  depends  on  the  phase  region  as  described 
in  the  last  section.  The  formula  for  CVjm  in 
a  mixed  phase  is  given  in  Appendix  II. 


Equation  of  state: 


n+1  n+1 

Pj+%  =  cp(vI>T>]j+% 


(5.13) 


where  p(vpT)  is  given  by  Eq.  (4.5). 

The  boundary  between  a  single  and  a  two-phase 
region  is  distinguished  by  the  transition 
pressure  p^,  at  which  the  relative  volume  of 
the  first  phase  is  given  by  v^. 
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5.2.  FORTRAN  Program 

The  FORTRAN  IV  program  used  for  the  calculations  is  given 
in  Appendix  III.  This  program  is  set  up  specifically  for  the 
problem  of  a  constant  driving  pressure  and  a  semi- infinite  slab 
of  material  so  it  does  not  treat  an  unloading  rarefaction  prob¬ 
lem  nor  the  reflection  of  a  shock  wave. 

The  heart  of  the  program,  insofar  as  we  are  concerned, 
is  in  the  calculation  of  a,  p  and  T,  following  that  of  v  and  q 
(Eqs.  (5.10)  through  (5.13)  above);  therefore  we  consider  it  in 
detail. 

When  the  material  of  cell  j  is  in  the  region  of  phase  1, 
as  determined  by  a  test  on  p,  the  numerical  computation  proceeds 
to  Eqs.  (5.12)  and  (5.13)  from  Eq.  (5.8).  When  the  pressure  in 
cell  i  reaches  the  transition  value  p^  and  overshoots  it,  tem¬ 
perature  and  the  relative  specific  volume  of  phase  1  must  be  re¬ 
adjusted  by  use  of  Eqs.  (5.11)  and  (5.12).  At  this  time  the 
special  program  parameter  NSA(J)  is  set  equal  to  2  for  the  cell 
j  and  kept  at  that  value  until  the  material  is  all  converted  to 
phase  2.  Calculation  of  pressure  is  accomplished  by  substituting 
v^  and  T  from  Eqs.  (5.11)  and  (5.12)  into  (5.13).  Once  the 
given  cell  becomes  pure  phase  2,  NSA(J)  is  set  equal  to  3,  and 
the  calculation  again  proceeds  as  for  phase  1.  However,  as  long 
as  we  use  the  relaxation  relation  (Eq.  (3.26)),  with  the  approx¬ 
imation  (v2_v^  =  constant),  we  always  have: 

(X  <  1 


v^  =  v  -  cy(v2  -  v^) . 


(5.14) 


If  we  are  reminded  that  pressure  is  essentially  calculated  from 
v^,  then  from  Eq.  (5.14),  regardless  of  whether  we  set  a  to 
unity  when  it  is  close  to  unity,  the  value  of  is  practically 
the  same.  Therefore  we  can,  for  the  present  case,  discard  the 
case  for  NSA(J)  =  3,  which  is  enclosed  by  the  dotted  lines  in 
the  flow  chart. 

The  expanded  flow  chart  for  this  part  of  the  program  is 
given  in  Fig.  5.3. 

5.3.  Numerical  Results 

5.31  Transient  Case 

The  particular  concern  here  is  with  development  of  the 
double  wave  structure  associated  with  the  phase  transition.  A 
uniform  pressure  is  applied  at  time  t  =  0  to  the  surface  of  a 
half  space,  x  =  0,  and  maintained  constant  as  the  plane  wavefront 
develops  and  propagates  inward.  For  an  applied  pressure  of 
0.200  Mbar*,  Figs  5.4  and  5.5  show  pressure  profiles  of  the  de¬ 
veloping  wave  at  times  measured  from  the  first  application  of 
pressure.  Fig.  5.4  illustrates  the  development  and  decay  of  the 
first  wave  caused  by  the  a  -  e  phase  transition.  After  about 
ten  relaxation  times  the  profile  has  the  clear  double-wave  struc¬ 
ture  shown  in  Fig.  5.5. 

The  thickness  of  the  first  wave  front  is  determined  by 
q  and  by  the  choice  of  space  increment,  Ax  (Fig.  5.1).  The  width 
of  the  second  front  is  controlled  by  the  relaxation  time  of  the 
phase  transition.  Nhen  it  is  well  separated  from  the  first  wave, 

*1  Megabar  =  10^  dynes/cm^ 
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Computation 

of 

U,  Q,  V 


Subroutine  2MIX 
Mixed  phase 
calculation  of 
p  and  T . 
(5.10),  (5.11), 
(5.12),  (5.13). 


i - 


l - 1 

I  r- — i  !  1 

|  -  NSA=3I — 7 

|  I  Test  —  I  1 

i  !  L _ j 

1  feil - 
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i _ 


i-' 


Set  NSA=*2 


Subroutine 
PHASE  1 
Single  phase 
calculation 
of  (p,T)  by 
(5.12  and 
(5.13). 


Test  (p,T) 

■1 

Crossed 

the 

1 

boundary 

{yesH 


|  j  i 

i - *1  Second  phase  subroutine 

•  i 


,L 


End 

of 

Cycle 


.  200  Mbar 
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the  relative  magnitude  of  q  in  the  second  front  is  negligible, 
as  shown  in  Fig.  5.6.  The  coefficients  C  are  varied  to  deter- 

la 

mine  effects  of  q  on  the  transient  wave  profiles. 

Richtmyer  (36)  has  shown,  using  quadratic  q,  that  a  small 

coefficient  of  pseudo-viscosity  produces  an  ocillatory  output 

even  when  the  stability  condition  (Ax/At)  is  satisfied.  Figs.  5. 

through  5.9  show  similar  changes  in  profiles  for  various  C^. 

Figs.  5.10  and  5.11  give  the  profile  at  fixed  times  for  three 

different  values  of  C^.  From  these  it  is  quite  clear  that  the 

shock  profile  converges  to  the  same  form  after  about  three 

relaxation  times.  Details  of  the  relaxation  process  at  early 

times  can  be  obscured  by  oscillations  when  is  too  small  as 

these  figures  show.  A  very  large  produces  so  much  damping 

that  sudden  changes  in  profile  are  prevented.  This  can  allow  q 

to  control  the  profile  of  the  second  shock  as  well  as  the  first. 

A  value  of  C  =0.1  was  found  satisfactory  for  most  of  the  cal- 
Li 

culations  described  here. 

Novikov  observed  a  drop  in  pressure  behind  the  first 
shock.  These  calculations  sometimes  show  such  a  drop,  but  it 
is  more  likely  due  to  oscillations  in  the  output  than  to  a 
physical  effect. 


.200  Mb 


oastf  C/i 


(ea'eqoix^)  i  ‘ejcnssaad 


Fig. 5. 8. — Wave  Propagation  v/ith  C  =  .1. 
md  At* are  space  and  time  increments  used  in 
numerical  integration. 
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md  are  space  and  time  increments  used 
:ne  numerical  integration. 
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Both  pressure  and  or  are  shown  in  Fig.  5.12  as  functions 
of  time  at  .01  cm  from  the  surface  for  3  values  of  CL*  The  be¬ 
havior  of  is  sensitive  to  C^;  it  reflects  the  fluctuations 
in  p.  or,  on  the  other  hand,  is  relatively  independent  of  C^; 
its  behavior  is  controlled  by  the  relaxation  time  t.  This 
suggests  that  the  decay  in  amplitude  of  the  first  wave  is  essen¬ 
tially  independent  of  the  artificial  viscosity. 

The  decay  of  the  precursor  is  shown  in  Fig.  5.13.  The 
rate  of  decay  at  early  times  is  closely  related  to  a  simple 
exponential,  as  shown. 

This  behavior  is  plausible  on  the  basis  of  the  follow¬ 
ing  model.  Eq.  (3.23)  can  be  written 

O 

dp/dt  =»  a  dp/dt  +  (m-jn2  ”  m2nP*^vlT ,Qf^G 

This  is  the  analog  of  Eq.  (9)  in  reference  (50).  By  the  same 
arguments  used  there,  we  can  arrive  at  the  analog  of  Eq.  (19) 
of  ref.  (50) : 

dp-^dt  «  -  (mln2  "  m2nl^aeq  “  °0/tG 

providing  the  path  of  the  number  one  shock  lies  along  a 
characteristic,  which  is  nearly  true.  Here  p-^  is  pressure  at 
the  peak  of  the  first  shock,  assumed  to  be  a  discontinuity, 
hence  it  lies  on  the  metastable  surface  v^(p,T).  With  this 
condition,  a  E  0.  Now  with  the  sweeping  approximations  that 
Vj-V£  -  Av  *  constant,  that  the  entire  process  is  temperature 
independent,  and  *hat  =  Cp2>  we  obtain  for  the  decay 
equation: 

dp^/dt  -  (&v  «  /2T)dp/dv^  . 


200  Mb 


Fig.5.12r-<1  and  Speed  of  Transformation  at  Fifth  Cell 


O  Precursor  amplitudes 


Fig.  5.13. — Decay  of  First  Shock  in  Iron  Resulting 
from  Fnase  Transition. 
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Here  vm  and  pm  are  pressure  and  volume  where  the  Hugoniot  first 
enters  the  mixed  phase.  These  expressions  yield 


P1  =  Pd  +  (x  L  v/2Ut)  dp/dvj^  ,  £  vra-Av 


3  pm  +  ^pD"Fm^  exp(-x/2Ut)  , 

v  -Av  £  V-,  *  v_ 
m  -L  m 


where  x  =  Ut.  Figure  (5.17)  shows  that  in  the  present  case 
(pD  =•  driving  pressure  =  200  kbar) ,  the  central  formula  applies, 
therefore  we  should  expect  to  find  that 


p-^  -  130  =  70  exp  (-x/2Ut) 

The  difference  between  this  curve  and  the  numerical  results, 
shown  in  Fig.  5.13,  is  due  to  non-linear  effects. 

In  Figs.  5.5  and  5.11  there  are  arrows  labelled  A  and  B. 
These  indicate  the  shock  front  position  which  would  be  predicted 
at  the  indicated  times  from  the  Rankine-Hugoniot  jump  conditions 
A  for  the  first  shock,  B  for  the  second.  The  difference  be¬ 
tween  this  predicted  arrival  time  and  the  one  obtained  in  the 


I 
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numerical  integration  arises  because  the  shock  profile  has  not 
reached  its  steady  state  or  permanent  regime.  For  the  permanent 
regime,  to  which  -jump  conditions  apply,  the  locus  of  p+q  vs  v 
lies  along  two  straight  lines,  OM  and  MC  in  Fig.  5.14.  M  repre¬ 
sents  the  break  in  pressure  between  the  first  wave  and  the 
second,  and  the  velocities  of  the  first  and  second  waves  are 
proportional  to  the  square  roots  of  the  slopes  of  these  two 
lines.  'Jhen  the  pressure  is  applied  suddenly  and  the  relaxation 
time  is  long,  the  state  of  an  element  will  rise  to  a  point  on 
the  extended  metastable  curve  of  phase  1,  say  C1  in  Fig.  5.14, 
and  then  proceed  toward  the  equilibrium  point  0  along  an  isobar. 
During  the  process  of  precursor  decay,  the  locus  will  lie  along 
intermediate  curves  such  as  OabC.  In  fact  then  the  leading  wave 
and  the  developing  profile  behind  it  propagate  at  cirst  with  a 
velocity  proportional  to  the  square  root  of  the  slope  of  the 
chord  OC1,  only  gradually  approaching  the  smaller  steady  state 
velocities.  Consequently  the  wave  front  position  calculated  by 
integration  of  the  flow  equations  should  always  lie  ahead  of  the 
position  predicted  by  the  jump  conditions  whenever  the  constitu¬ 
tive  relations  include  rate  or  time  dependent  processes.  Loci 
for  pfq  obtained  in  the  calculations  for  several  positions  and 
two  values  of  are  shown  in  Fig.  5.15.  They  do  indeed  dis¬ 
play  the  behavior  described  in  Fig.  5.14;  moreover  the  locus  is 
relatively  independent  of  C^. 

It  is  possible  to  calculate  the  equilibrium  Hugoniot 

curve  for  a  material  undergoing  a  phase  change  using  only  the 

an 

jump  conditions  andAequation  of  state  (17).  However,  once  a 


Fressure,  (p+q) (Kilobars) 
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program  is  available  for  integrating  the  flow  equations,  as  in 
this  case,  it  turns  out  to  be  easier  to  run  the  program  with  a 
uniform  driving  pressure  and  tabulate  p  and  v  in  the  uniform 
region  far  behind  the  shock  fronts  than  to  do  the  equilibrium 
computation.  That  has  been  done  and  the  results  are  shown  in 
Fig.  5. 16. 

TI  is  the  temperature  independent  solution  (which  is  on 
the  second  phase  isotherm) ,  and  TDl  and  TD2  are  temperature  de¬ 
pendent  solutions  with  different  isothermal  compressibilities 
(a£  -  3.4  and  a2  =*  2.4,  respectively).  The  difference  between 
TI  and  TDl  is  due  only  to  the  temperature  dependence  in  the 
equation  of  state.  The  reason  TI  lies  above  TDl  can  be  explained 
as  follows:  since  dp/dT  is  negative  in  the  coexistence  region, 
the  transition  produces  a  temperature  decrease.  This  decrease 
is  found  to  be  larger,  except  near  the  transition  point,  than 
the  temperature  rise  in  the  first  shock  (about  20°) ,  hence  it 
gives  a  slight  pressure  drop  for  the  temperature-dependent 
equation  of  state.  Experimental  values  measured  by  Minshall  (13) 
are  indicated  by  crosses.  The  differences  between  these  and 
the  calculated  curves  are  substantial.  It  is  quite  possible 
that  a  curve  passing  through  points  B  and  C  can  be  developed  by 
allowing  v^-v^  to  vary  with  p.  Point  A,  however,  appears  to  be 
unattainable  within  what  are  here  believed  to  be  reasonable 
limits  of  the  thermodynamic  and  transition  parameters. 


Pressure,  p  (Kilobars) 
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Fig.  5.16. — The  Hugoniot  Curve  in  and  beyond 
the  Coexistence  Region. 
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5.32  Permanent  Regime 

After  the  shock  of  the  last  section  has  propagated  far 
into  the  medium,  the  profiles  of  the  first  and  second  shocks 
are  expected  to  become  unchanging,  though  they  continue  to  sep¬ 
arate  because  of  differing  velocities.  When  this  happens,  these 
profiles  should  be  desciibed  by  a  permanent  regime  solution  to 
the  flow  equations  (37).  Such  a  solution  is  obtained  here  in 
order  to  determine  how  far  the  wave  must  travel  to  closely 
approach  the  permanent  regime  and  to  provide  an  independent 
check  on  the  numerical  integration  for  the  transient  case.  We 
proceed  by  setting  (*/*t)x  =  0  for  all  variables  in  Eqs.  (3.3)- 
(3.5)  and  (3.43).  The  resulting  equations  are,  for  the  temper¬ 
ature  independent  case: 


pu  du/dx 

=  -dp/dx 

(5.15) 

d(pu)/dx 

ii 

o 

w  • 

*o 

c 

II 

3 

(5.16) 

udcv/dx 

=  (ffeq-<*)/  T 

(5,17) 

v  =  v^ 

+  (v2-v1)q- 

(5.18) 

2(P)  "  vi 

(p)  -  const 

(5.19) 

p  = 

P(vx) 

(5.20) 

Combining  Eqs.  (5.16)  and  (5.15)  yields  the  Earnshaw  relation: 

p-p  =  PV  (v  -v)  =  m^(v  -v)  (5.21) 

O  0  0  o 

where  U  is  shock  velocity.  Combining  Eqs.  (5.17),  (5.18),  (5.19) 
and  (5.21)  yields  an  equation  for  dp/dx  in  the  transition  region 


dp/dx 


m  (v^*^) (o’e^-o')/T(l-hn^dv^/dp)v  . 


(5.22) 
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From  the  definition  of  »et*  we  have 

Q,eq  ■  0>  p  s  pQ 

*  (v“v0)/ (va“v0)  *  po  *  p  *  pa  (5*23^ 

*  1  >  Pa  s  P 

where  pQ,  pa,  v  ,  vQ  are  defined  in  Fig.  5.18.  Combining 
Eq.  (5.23)  with  (5.18),  (5.21),  and  (5.22)  yields  the  following 
expressions  for  dp/dx: 

dp/dx  -  m3(v0-v1)/T(l-hn2dv1/dp)(m2v0-p+po)  , 

po  S  p  *  pa  (5.24) 

dp/dx  =*  m(m2(2v0-v-L-va)  +  po-p)/T(l+m2dv^/dp) 

•  (m2v0-p+po)  ,  pa  s  p  s  pf  (5.25) 

Here  is  a  known  function  of  p,  Eq.  (5.20).  Examination  of 
Eqs.  (5.24)  and  (5.25)  shows  that  dp/dx  =  0  at  p  =  pQ  and 
p  =  p^.  In  Eq.  (5.22),  m  <  0  for  a  forward- facing  shock, 
tyec*  -  or  >  0,  v^-V2  >  0,  and  1  +  m2dv^/dp  >  0,  from  Fig.  5.18. 
Therefore  dp/dx  <  0  throughout  the  transition.  The  qualitative 
features  of  dp/dx  are  shown  in  Fig.  5.18. 

Eqs.  (5.24)  and  (5.25)  have  been  integrated  for  the 
iron  transition,  assuming  that  no  temperature  changes  occur  in 
the  shock.  The  results  are  compared  with  the  temperature- 
independent  transient  case  in  Fig.  5.19  for  a  driving  pressure 
of  .200  Mbar.  The  shock  width,  defined  as  Ax  * 4p/ |ldp/dx|raax, 
is  for  this  case  about  .3  cm.  The  steady  velocity  of  the 
second  shock  with  respect  to  the  material  ahead  of  it  is  Ujj  * 
0.37  cra/psec.  The  relaxation  time  assumed  for  the  calculation 
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Fig.  5.17.— Path  of  the  Permanent  Regime  Solution 


Fig. 5. 18r-dp/dx  for  Permanent  Regime  Solution 


Fig.  5. 19. -“Comparison  of 
permanent  regime  with  transient 
profile  of  second  shock. 

Q  Run  #67108-2,  Ax  *  .01, 

At  -  .0025  Usee,  t  -  15t 

o  Run  #67108-3,  Ax  »  .001, 

At  -  .0005  lisec,  t  -  2.25 

—  Permanent  regime  profile 


Fig. 5. 19.  — Comparison  of  Permanent  Regime  with  Transient 
Profile  of  Second  Shock. 
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is  t  ■  1/3  nsec,  so  Ax  «  2.4  U^^t.  Profiles  obtained  from 
the  transient  calculation  at  t  -  2.25t,  15t,  22. 5t  are  shown 
for  comparison.  The  agreement  is  very  good  at  15t.  This 
agreement  is  additional  evidence  of  the  validity  of  the  inte¬ 
gration  procedure  for  the  transient  case. 
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APPENDIX  I 


HEAT  CAPACITY  OF  MIXTURE 
Equation  (31)  of  Reference  (7)  gives 

CVM  "  “  T<*V/&P>SM  (dp/dT>2  (D 

where  (Sv/dP)^  is  the  slope  of  an  equilibrium  adiabat  in  the 
mixed  phase  region.  Equation  (16)  of  the  same  reference  gives 

(5V/ap)sM  *  dVdp  “  (dSx/dT)(dT/dP)2 

-  (V  -  VjHdT/dP)2  d2P/ dT2  .  (2) 

Using  the  identities 

T  dS^/ dT  =*  Cyl  +  T(bP/5T)yl  dV^dT 

and 

dV1/ dT  =*  (bVj^/d^p  +  (av1/3P)T  dP/dT  , 

Equations  (1)  and  (2)  can  be  transformed  to  yield 

-  CV1  -  T(bP/av1)T  (dV1/dT)2  +  T(V-V1)  d2P/dT2  (3) 

The  first  two  terms  of  Equation  (3)  correspond  to  Equation  (85) 
of  Reference  (38).  The  third  term  may  be  important  if  the  state 
point  is  not  near  the  boundary  of  phase  I. 
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APPENDIX  II 

PROGRAM  FOR  CALCULATING  WAVE  FLOW  IN  ONE  SPACE 
DIMENSION  FOR  PLANE,  CYLINDRICAL,  OR 
SPHERICAL  GEOMETRY 

A.  Program  Name:  BURN 

B .  Program  Description 

This  program  is  an  adaptation  of  one  written  by  John  0. 
Erkman  at  Stanford  Research  Institute.  It  integrates  the  equa¬ 
tions  of  flow  through  one  space  and  one  time  dimension  for  ar¬ 
bitrary  initial  and  boundary  conditions.  Integration  is  carried 
through  shock  fronts  by  means  of  an  artificial  viscosity. 

Initial  and  boundary  conditions,  input  parameters  and  output 
statements  are  contained  in  subroutines  so  they  may  be  readily 
altered. 


Subroutine  DECIDE  contains  the  input  parameters  which 
define  the  problem.  Comments  in  the  Listing  (Section  C)  should 
make  this  subroutine  self-explanatory,  except  perhaps  for  the 
following: 

S  is  an  integer  index  used  to  label  the  various  regions 
in  the  problem,  S  =  2, SI 


BUKN(S)  is  an  index  which  defines  the  material  in  region  S: 
for  example  if  the  geometry  Is  spherical  and  the  problem 
consists  of  a  sphere  of  explosive  surrounded  by  an  A1 
shell,  SI  *  3,  BURN (2)  =  1,  BURN(3)  -  4.  Other  values 
of  BURN(S)  are  defined  in  the  program  listing. 


0PTI0N  is  an  index  defining  the  driving  system  for  the 
problem.  0PTI0N  *  1,2,3,  or  4  for  a  pressure  pulse 
applied  to  the  left  boundary. 
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TAU  is  a  characteristic  time  parameter  for  the  applied 
pressure.  For  its  exact  meaning,  consult  the  listing 
of  the  MAIN  program.  It's  units  are  microseconds. 

LEFTP  is  in  Megabars  (Mb) 

TQUIT,  microseconds 

Values  of  Z0N(S)  and  L(S)  need  be  given  for  S  =  2,  SI. 

The  primary  output  consists  of  tables  of  values  of  particle 
velocity,  pressure,  etc.,  vs  J  for  each  time  and  cycle 
indicated  in  DESCRIBE. 


C.  Program  Listing  and  Sample  Output 


ooo  ooor>  o  or>oonr>o 
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//BOON  6XEC  FORTGCLG 
//FOPT.SYSIN  00  * 

THIS  IS  A  3NE-0 IMFNSI ONAL  Q-CODE  ADAPTEO  FROM  ONE  WRITTEN  BY 
JOHN  0.  EPKMAN  OF  SRI,  WHO  MODFLED  AFTER  ONF  WRITTEN  BY 
MARX  WILKINS  OF  LRL .  AN  APPROXIMATE  DESCRIPTION  CAN  BF  FOUND 
IN  "COMPUTATIONAL  PHYSICS.  VOL.  Ill,"  EDITED  BY  ALDER  AND 
FFRNBACH  ANO  ROTENBERG 

THE  PROBLEM  IS  SET  UP  IN  SUBROUTINE  "DECIDE." 

THE  LEANING  OF  KF Y  SYMBOLS  IS  DESCRIBED  THFRE . 

COMMON  /C17.0N/  H(9),BURN(9),L(91,DX(9) ,Sl,RM0<9) 

COMMON  /C2TIME/  T  IMES ,C YCLE , DEL T , DTN,D TMX, TL I MA( 300 ) , JCRI T, 

1  TQU I T, TAU 

COMMON  /C3CTRL/  COUNT S, JSTAR ,JPE , JPR,JQUI T, LAST, CYCLES 
COMMON  /C4FL0W/  J ( 300  ) , VI  300  I , X( 300 1 , 0 ( 300) , P ( 300  I , E ( 300 ) , OA , VN , 
1MASSI 300) ,CSPI 300 ) 

COMMON  /C7GNRL/  ALP , OPT  I  ON, CONA , CO , LEF TP 

INTFGER  H, BURN, S, SI, ZON, CYCLE, COUNTS, CYCLES, ALP, OPT  ION, H2,HS1»HS* 
1  BURNS, HS2 

REAL  L, MASS, LINEAR, LEFTP 
CALL  DECIDE 

THE  FOLLOWING  00  LOOPS  ENDING  AT  STATEMENT  9  CALCULATE  THE 
POSITION  OF  THE  J'TH  CELL  IN  CM  AND  ITS  MASS  IN  GRAMS.  RHO(S)* 
DENSITY  OF  RFGION  S  IN  GRAMS/CC. 

DO  9  S=2 , SI 
HS1  =  H( S-  1  )  +  l 
HS2=H( S) 

DO  9  J=HS1 , HS2 
X( J+l ) =X( J) «-DX( S) 

Q  MASS!  J)  =  (X(  J  +  l )**ALP-X( J)**ALP)*RHO(S> 

THF  VARIABLES  IN  THE  FOLLOWING  F3UR  WRITE  STATFMENTS  HAVE  BEEN 
DEFINED  IN  SUBROUTINE  DECIDE. 

WRITE (6, 951 » ALP, OELT,OTMX, CONA, CQ 
WRITE (6, 952 J C YCLF S, COUNTS , JQUI T 

9  52  FORMAT! ’O' , ’CYCLES’ ,6X, 'COUNTS ' , 6X , » JQUI T ' / 1  6 ,4X , I  6, 6X , ,  U 
WRITE! 6, 957) SI , (BURN! S) ,S=2,Sl) 

WRITE (6, 961 ) TAU, LEFTP, U!1 ), OPTION 
C 

IF  (OPTION. EO. 6)  CALL  FLIER 
IF  (OPTION. NE. 6)  JSTAR=5 
C 

CALL  WRITE l 
C 

CqSQ=CQ**2 
CQSQ4=4.0*CQSQ 
L  INEAR  =  1 .O+CONA+CONA 

C  "TIMES"=T,  THE  TIME  VARIABLE,  MFASUREO  FROM  ZERO. 

T IMES=0. 0 
CYCLE=0 

C  " JCR l T"=NO.  OF  SPACE  CELL  FOR  WHICH  TL IMA! J)  HAD  ITS  MINIMUM 

C  VALUE  IN  THE  PREVIOUS  CYCLE. 

JCR IT=0 


oooooooo  ooo 
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C 

r 

C 


C- 

r. 


40 


so 


SI 

s? 

S3 

60 


70 

78 

79 
R0 


THE 

AFTER  THE 


IS  AN  INTEGER  CONTROL  PARAMETER  which  DIRECTS 

MWRITE"  T0  PERFORM  a  terminal  computation 
INTEGRATION  HAS  BEEN  COMPLETED. 

L AST=0 

0TN*DELT  THE  VALUE  °F  "°ELT"  CALCU*-ATE0  IN  THE  CYCLE  BEFORE  LAST. 

DEL  T l  -DEL7*-DEL  T 
-START  OF  TIME  LOOP 

PPEAK=0?0IHUM  VALUE  °F  PRE$SURE  CALCULATED  IN  PREVIOUS  CYCLE. 

TLIMB=TLIMA( JCRIT JsMINIMUH  »ALUE  OF  TLIMA(J). 

TL I  MB3 1.0 

times=times+oelt 
cm  f=cycle*i 

J=1 

s«? 

J  1=2 
JT=3 

PLFFT=0.0 

-r?nNI«;fV^UcTF  P  F°R  FIRST  CELL  AN0  U  Am  *  ™  LEFT  BOUNDARY 

oO  TO  (51,52,53,54,60,60) .OPTION 
(FITIMES  .LE.  TAU)  PLEFT=LEFTP 
GO  TO  60 

IFITIMES  .LE.  TAU)  PLEFT-(  (-TIMES /TAU)  *1.0) *LFFTP 
GO  TO  60 

PLFFT«LEFTP*FXP(-0.46*TIHES ) 

CONTINUE 

•START  OF  J-LOOP 
IF(  J.GT.H(S)  )  S  =  S  *1 

DENU=(X( JT)-X( Jl) )/V(Jl)*(X(JI)-X(J))/V(J) 

U(  J1 )  MOELT  I*(  P(  J  )-P(  Ji)«.Q(  j  j_q(  ji)  j  j  /qenu*U(  Jl ) 

X  ( J I  =  XA 

XA=DELT*U(JI)*X(J1) 

IF(J  »  EQ.  H(SD)  X  (  J  l  )=XA 
I Ft ABS*U( J 1 ) )  .LT.  5.0E-5)  U(J1)=0.0 
VN=(XA**ALP-X( J)*»ALP)/MASS(J) 

DELU=U(  J  l )  — U (  J  ) 

OELX=XA-X( J  » 
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-GET  0  FOR  SHOCK 

•*QA”=NEW  VALUE  OF  ARTIFICIAL  VISCOUS  STRESS;  Q(J>  IS  "OLD'*  VALUE 

QA3-DELU*(CQSQPA8S( OELU)+CONA*CSP(J) )/VN  L  VALUE‘ 

IF(QA  ,LT.  0.0)  QA=0.0 

MTL  IMA  (  J  l*T  twp  P4BAMCTCO  uccn  >u  /••.»•...  ......  _  . . .  . 

fnn  rue  MCYT  True  ere  A""  *’  '"«*■*•«  l- «  '  « "H»  I  Mt  VALUE  OF  DELT 

FOR  THE  NEXT  TIME  STEP.  AT  THIS  POINT  VN,  THE  NEW  VAIitF  df 

»mSr5ufE 0M  THIS  ”  ,s  ,‘UNSFI;RRE0  ™  ™e 

TL  FMA(  J)=0ELX/(LINEAR'i'CSP(J)*CQSQ4*ABS(DELU)  ) 

CALL  EQSNSrJ) 


o<-oo  <_  o  u  o 
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I F ( ABS( P C  J I )  .LE.  1.0E-5)  P( J! =1 . 01 34E-6 
VI  J)=*VN 
0(  J  )  =  QA 

IF(TLIMA(J).GE.TL!MB)  GO  TO  129 
J  CR I T=J 
TLlMB=TLIMA(  J) 

1 2P  IF(P(J)«-Q(J).LE.PPEAK>  GO  TO  131 
PPFAK=(P(J)*Q(J)) 

C  "JPMAX"  IS  THE  VALUE  OF  J  FOR  WHICH  PPEAK  TAKFS  ITS  MAXIMUM  VALUE. 
JPMAX=J 
131  J=J*1 
Jl=J+l 
UT=Jl+l 

" JST AR"  IS  A  VALUE  OF  J  AHEAD  OF  THF  WAVE  FRONT  AT  WHICH  COMPUTING 
STOPS  PREPARATORY  TO  ADVANCING  TIME  BY  ANOTHER  INCREMENT.  JSTAR 
IS  ADVANCED  BY  UNITY  WHENEVER  PARTICLE  VELOCITY,  UIJSTAR? , 

BECOMES  NON-NEGL I GI 8LE. 

I F ( J.LE.JSTARU)  GO  TO  70 

- TEST  TO  ADVANCE  JSTAR 

IF(ABS(U(JSTARni  ).GT.1.0E-5»  JSTAR=JSTAR*1 

lF((CYCLE.EQ.CYCLf  1. OR. (TIMES. CE.TQUIT1 .OR. { J. EQ. JOUI T) 1  GO  TO  169 
IF(CYCLE.GT.10)C0  .  20 

TO  CHANGE  FREQUENCY  UF  Pi  NT-OUT,  A  STATEMENT  CAN  BE  INSERTED 
HERE:  "IF  (CYCLE. GT.  K)  COUNTS=MN"  WHERE  "K"  AND  "MN"  ARE 

INTEGERS  TO  BE  CHOSEN  BY  THE  PROGRAMMER. 

IF(MDD(CYCLE, COUNTS!. NE.O)  GO  TO  180 
GO  TO  170 

169  L  AST= 1 

170  JPB= 1 
JPF=JSTAR+2 
CALL  WRITE 

1  BO  DTNH1=0.6*TLIHB 

IE? DTNH1/DELT .GT. 1. 1)  DTNH1=1 . 1 *DELT 
IF(0TNH1.GT.DTMX»  DTNH1=0TMX 
OTN=DELT 
DELT^DTNHl 
DELT I=DTN*DFLT 
GO  TO  40 

9«51  F0RMAT(1H1,6X,3HALP,9X,4HDELT,  1 1 X, 4HDTMX,  11 X , 4HC0NA ,  1 2X , 2HCQ/7X ,  1 1 
-, 4F15.6) 

9S7  FORMAT(lHO,8X,2HS1.5X,8HBURNU..m»I2,9X, 91 5/2X1 
9  61  FORMAT ( 1H0, 5X,  '.HT AU,  1 6X,  5HLEFTP,  14X,4HU(  1 ) ,  15 X, 6H0PTI 0N/3E19. 8, 1 8 ) 
FND 

SUBROUTINE  DECIDE 

COMMON  /C1Z0N/  H( 9) ,BURN( 9) ,Ll 9 ) , DX( 9) , SI ,RMO (9) 

COMMON  /C2TIMEZ  T IMES , CYCLE, DELT, DTN.DTHX, TLI MA( 300) , JCRI T , 

1  7QU IT, TAU 

COMMON  /C3CTRL/  COUNTS,  JSTAR ,JPE , JPB,JQUI T, LAST , CYCLES 
COMMON  /C4FL0W/  U( 300 ) , V( 300! , XI 300) ,Q( 300 1 , P (300 ! , E ( 300) , QA, VN, 

1  MASS (300 ) ,CSP ( 300 1 

COMMON  /C7GNRL/'  ALP, OPTION, CONA,CQ,LEFTP 
DIMENSION  ZON ( 9 ) 


oooooooooooooooo  ri  o  o  o  o  o  (  ^  ooooooooor*»o  o  nnnooo^onooo 
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INTEGER  H,BURN,S,S1, ZON, CYCLE, COUNTS, CYCLES, ALP, OPT  ION, H2.HSl.HS, 
1  BURNS, HS2 

REAL  L, HASS, LINEAR, LEFTP 


CHOOSE  GEOMETRY.  ALP  IS  AN  l NTEGER  LABEL  WHICH  IS  TO  BE  SET 
ACCORDING  TO  THE  GEOMETRY  OF  THE  PROBLEM. 

X  IS  THE  EULER l AN  SPACE  COORDINATE.  THE  INITIAL  VALUE  AT  T=0 
OF  THE  LEFT  BOUNDARY  OP  CELL  1  IS  SET  HERE.  THE  POSITIONS  OF 
OTHER  CELL  BOUNDARIES  ARE  CALCULATED  IN  MAIN  FROM  THE  NUMBER  OF 
ZONES  AND  THE  DIMENSIONS  OF  THE  PROBLEM. 

FOR  CYLINDRICAL  AND  SPHERICAL  PROBLEMS,  LFFT  BOUNDARY  IS 
INTERPRETED  AS  INNER  BOUNDARY. 

ALP  =  1  FOR  PLANE  GEOMETRY 
ALP  =  2  FOR  CYLENDRICAL  GEOMETRY 
ALP  =  3  FOR  SPHERICAL  GEOMETRY 


ALP=  1 

- CHOOSE  COORDINATES  OF  FIRST  CELL 

X(1)=0.0 

- NUMBER  OF  REGIONS  PLUS  ONE  (NOT  TO  EXCEFD  9) 

THIS  PROGRAM  CAN  BE  RUN  WITH  SEVERAL  REGIONS  OF  DIFFERENT 
MATFRIALS.  THE  NUMBER  OF  SUCH  REGIONS  IS  DENOTED  BY  AN  INTEGER 
Sl-I.  THIS  PECULIAR  CONVENTION  ARISES  BECAUSE  OF  A  CHARACTERISTIC 
OF  FORTRAN— ZERO  INDICES  ARE  NOT  ALLOWED.  EACH  DISTINCT  REGION  IS 
DENOTED  BY  AN  INTEGER  LABEL  S.  S=2  IS  THE  LEFT-MOST  REGION,  S=3 
THE  NEXT  TO  THE  RIGHT,  ETC.  UP  TO  SI. 

EACH  REGION  IS  DIVIDED  INTO  A  NUMBER  OF  SPACE  ZONES  OR  CELLS, 
ZONISI,  THE  NUMBER  OF  CELLS  UP  TO  AND  INCLUDING  REGION  S 
(STARTING  WITH  THE  LEFTHOST  REGION)  IS  H( S)=SUM( ZON( K) ) , 

K=2  TO  S,  INCLUSIVE. 

S  1=2 

- MATERIAL  IN  REGIONS 

"BURN(S)"  IS  AN  INTEGER  LABEL  WHICH  DEFINES  THE  MATERIAL  OF 
REGION  S. 

BURN (S)  =  1  FOR  EXPLOSIVE 
RURN(S)  =  2  FOR  VOID 
BURN ( S )  =  3  FOR  LIQUID 
BURNIS)  =  4  FOR  SOLID 
BURN(S)  =  5  FOR  PHASE  TRANSITION 
BURN { 2 )-5 
- SET  OPTION 

"OPTION"  IS  AN  INTEGER  LABEL  WHICH  DFSCRIBES  THE  TYPE  OF  PROBLFM 
TO  BE  SOLVED.  IF  OPTION=l,2,  OR  3,  THE  PROBLFM  IS  ONE  IN  WHICH 
A  SPECIFIED  PRESSURE  IS  APPLIED  TO  THE  LEFT  HAND  BOUNDARY.  IF 
OPTlON=5,  AN  EXPLOSIVE  REGION  IS  INCLUDED  AND  ITS  DETONATION 
PROVIDES  THE  DRIVING  FORCE.  0PT!0N=6  MEANS  THAT  THE  FIRST  REGION 
(  S*2 )  IS  A  FLYER  P' ATE  WHICH  HAS  JUST  COLLIDED  WITH  THE  SECOND 
REGION  ( S*3 )  AT  THE  START  OF  THE  PROBLEM.  WHEN  THIS  HAPPENS, 

EACH  CELL  IN  REGION  1  «S«2i  IS  GIVEN  THE  FLYER  PLATE  VELOCITY 


U(l),  EXCEPT  THE  ONE  ADJACENT  TO  REGION  2(S=3)  THIS  CFLL  AND  THfc 
FIRST  CELL  IN  REGION  2  <S*3)  ARE  GIVEN  VELOCITIES  U(l)/2  FOR 
SMOOTHING  PURPOSES.  WHEN  OPTION*!,  THE  TIME  DURATION,  TAU,  OF 
THE  APPLIFO  PRESSURE  MUST  BE  SET.  FOR  A  CONSTANT  PRESSURE  APPLIED 
AT  T =0,  SET  TAU  EQUAL  TO  A  LARGE  NUMBER,  SAY  500  (MICROSECONDS). 
FOR  0PTI0N*2,  TIME  TAU  IS  THE  TIME  AT  WHICH  THE  APPLIED  PRESSURE 
EQUALS  ZERO  IN  A  LINEAR  RAMP.  OPTIONS  HAS  A  BUILT-IN  TIME 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

f. 

C 

C 

c 

c 

c 
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c- 

c- 

c- 

c 

c 

c 

c 

c 

c 

r 

V 

c 

c 

c. 


0 
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CONSTANT.  THE  PEAK  APPLIED  PRESSURE  IN  EACH  CASE  IS  "LEFTP"  IN 
MEGABARS. 

OPTION  =  1  GIVES  SQUARE  PULSE 

OPTION  =  2  GIVES  LINEAR  PULSE 

OPTION  =  3  GIVES  EXPONENTIAL 

OPTION  =  4  UNASSIGNED 

OPTION  =  5  GIVES  NO  PULSE  IFOR  EXPLOSION! 

OPTION  =  6  GIVES  NO  PULSE  (FOR  FLYER  PLATE! 

PT I PN= 1 

IF  OPTION  =  1  OR  2»  SET  TAU 
T  AIJ=500. 0 

IF  OPTION  =  6 *  SET  UC 1 >  I  FOR  FLYER  PLATE! 

U(1)=0.0 

IF  OPTION  =  1*2,  OR  3,  SET  LEFTP  (PRESSURE  ON  LEFT  BOUNDARY! 

L  EFTP  =  0. 200 

VISCOSITY  COEFFICIENT  ICQ  FOR  QUADRATIC  AND  CONA  FOR  LINEAR! 

CON A=0. 1 
CQ=2«  0 

LENGTH  OF  RUN  MAY  BE  DETERMINED  BY  SETTING  ANY  OR  ALL  OF  NEXT 
WHEN  CYCLE=CYCLES  OR  TIMES=TQUIT  OR  J=JQUIT,  COMPUTATION  WILL 
STOP,  WHICH  EVER  OCCURS  FIRST. 

J  IS  THE  INTEGER  LABEL  OF  THE  SPACE  CELLS.  J=1  AT  THF  LEFTMOST 
CELL  OF  THE  LEFTMOST  REGION  AND  RUNS  TO  H(S1>,  THE  RIGHTMOST  CELL 
OF  THE  RIGHTMOST  REGION. 

CYCLES*  NUMBER  OF  INCREMENTS  IN  TIME 
TQUIT  (PROPAGATION  TIME) 

JQUIT  (NUMBER  OF  LAST  CELL) 

CYCL  FS=100 
TQUIT=?60 
JQUIT=25 
JQUIT=250 

THE  NUMBER  OF  ZONES  IN  REGION  K  IS  ZONIK) 

ZON ( 2 )  =50 

-THF  THICKNESS  OF  REGION  K  IN  CM.  IS  LIK) 

L (2»=5.Q 

-DELT  IS  STARTING  VALUE  FOR  DELTAT 

"DELTAT"  IS  THE  T IME- INCREMENT  FROM  ONE  CYCLE  TO  THE  NEXT, 
MICROSECONDS. 


1? 


DELT= .05 

-DTMX  IS  UPPER  LIMIT  FOR  DELTAT 
DTMX= • 05 

-PRINTOUTS  OF  CYCLES  IS  MODULO  COUNTS 

"COUNTS"  CONTROLS  PRINTING,  IF  C0UNTS=5,  THE  STANDARD  FLOW 
VARIABLES  U,P,Q,E,V,  ETC.  ARE  PRINTED  OUT  EVERY  FIFTH  CYCLE,  ETC. 

*i  IN  "DECIDE",  ThEN  AFTER  THE  FIRST 


10  oc  i 


FEW  CYCLES  INCREMENTED  TO  10  OR  20. 

COUNTS*! 

H(S)  IS  AN  INTEGER  LABEL  EQUAL  TO  THE  NUMBER  OF  SPACE  CELLS  TO 
THF  LEFT  OF  AND  INCLUDING  REGION  S. 

H  ( 1 )  *  0 
DO  12  S=2, SI 

DX(  SI  =  L ( S ) /FLOAT ( ZONI  S ) I 
H(SI*H(S-1I  ♦  ZON(S) 

-CALL  ROUTINES  TO  SET  INITIAL  REGIONS 

AT  THIS  POINT  CONTROL  IS  TRANSFERRED  TO  B_INIT(S)  FOR  S*2  TO  SI, 


n  o  o  o  o 
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WHERE  _  IS  AN  INTEGER  FRQH  1  TO  5,  CORRESPONDING  TO  THc  VALUE 
OF  BURN(S).  FOR  EXAMPLE  IF  Sl  =  4  AND  BURN(2I=1,  8URNI3)=4,  AND 
BURN(4I=3,  THEN  81INIT12),  B4INITI3),  B3INITI4)  ARE  CALLED  IN 
TURN?  I.E.,  THE  NEXT  THREE  STATEMENTS  WOULD  BE:  "CALL  BlINITI?)", 
"CALL  84  IN  IT  I  3 ) " ,  "CALL  B3INITI4)". 

CALL  85INITI2) 

RETURN 

102  FORMAT  (1014) 

906  FORMAT (2313 ) 

907  FORMAT (7I14F5.2/)  » 

END 

SUBROUTINE  EQSTIS.J) 

COMMON  /C1Z0N/  H(9),BURN(9),L(9),DX19) ,S1,RH0(9» 

COMMON  /C2TIME/  T IMES ,C YC LE , DELT, DTN.DTMX, TLI MA( 300 ) , JCRI T , 

1  TQUIT,  TAU 

COMMON  /C4FL0W/  U ( 3 00 ) , VI  300 1 , X I  300 ) , Q t 300) , P 1 300 ) , F I  300 ) , QA, VN, 

1  MASS! 300) »CSP( 300) 

INTCGER  H.BURN.S.Sl,  ZON,  CYCLE, COUNTS, CYCLES, ALP, OPT  ION, H?,HS1, MS, 
1  BURNS, HS2 

HURNS=8URN(S) 

GO  TO  (101,102,103, 104,105,106,107,108,109) , BURNS 

101  CALL  B1EQSTI S, J ) 

102  RETURN 

103  CALL  B3E0STI S, J  ) 

RETURN 

104  CALL  B4EQST I S , J ) 

RFTURN 

105  CALL  85EQST ( S, J ) 

RETURN 

106  RETURN 

107  RETURN 

108  RETURN 
104  return 

END 

SUBROUTINE  FLIER 

COMMON  /C1Z0N/  HI  9) , BURN! 9>,L(9) ,DX(9) »S1,RH0(9) 

COMMON  /C3CTRL/  COUNTS, JS TAR , JPE , JPB, JOU I T, L A  ST , C YCL E S 

COMMON  /C4FL0W/  U 1 300 ) , VI 300 ) , XI 300 ) ,3 1  300) , P 1 300) , E I  300) , OA, VN, 

1  MASS! 300) ,CSP(  300) 

INTEGER  H, BURN, S, SI, ZON, CYCLE. COUNTS. CYCI  FS. ALP.QPT !0N,H2; HS1 , HS , 
1  BURNS, HS2 

REAL  L, MASS, LINEAR, LEFTP 

J  STAR»H( 21+2 
H2=H( 2) 

00  43  J=1,H2 
43  U(J«-1)=»U(1) 

U(H2U}=0.5*UIH2*1) 


o  o  o  r>  n 
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RETURN 

END 

SUBROUTINE  B1INIT  <S) 

C  THIS  SUBROUTINE  IS  A  DUMMY  WHICH  ALLOWS  FOR  FUTURE  EXPANSION 
GO  TO  12 

ENTRY  B1E0ST(S,J) 

GO  TO  101 
1?  CONTINUE 
101  CONTINUE 
RETURN 


END 

SUBROUTINE  83INIT  (S) 

C  THIS  SUBROUTINE  IS  A  DUMMY  WHICH  ALLOWS  FOR  FUTURE  EXPANSION 
GO  TO  14 

FNTRY  B3  EOST ( S»  J  I 
GO  TO  121 
14  CONTINUE 
121  CONTINUE 
RETURN 
END 

SUBROUTINE  B4INIT  (SI 

C  THIS  SUBROUTINE  IS  A  DUMMY  WHICH  ALLOWS  FOR  FUTURE  EXPANSION 
GO  TO  13 

ENTRY  B4EQST ( St  J I 
GO  TO  90 
13  CONTINUE 
90  CONTINUE 
RETURN 
ENO 

SUBROUTINE  B5 1 N I T { S  ) 

THIS  SUBROUTINE  IS  WRITTEN  SPECIFICALLY  FOR  IRON  WITH  A  SHOCK- 
INDUCED  PHASE  TRANSITION. 

THE  PARAMETERS  ARE  DEFINED  IN  APPENDIX  II  OF  "EQUATION  OF  STATE 
IN  SOLIDS,"  BY  G.  E.  DUVALL,  G.  R.  FOWLES,  AND  Y.  HORIE,  SUMMARY 
REPORT  ON  CONTRACT  NO.  DA-04-200- AMC-1 702 ( X ) ,  BALLISTICS  RESEARCH 
LABORATORY,  A8ER0FFN  PROVING  GROUND,  MD. ,  FEB.,  1967. 


C 

C 


COMMON  /C1Z0N/  H(9),BURN(9)*L(9)»DX(9)?S1»RH0(9) 

COMMON  /C2TIME/  T IM ES , CYCLE , DEL T, DTN ,DTMX , TLI MA(300) , JCRIT, 
l  TOU  IT, TAU 

COMMON  /C4FL0W/  U (300  ) , Vt 300), X( 300) ,3(300) ,P (300) ,E( 300) ,QA,VN, 
1  MASS( 300) ,CSP( 300) 

COMMON  /C5THER/  TMP ( 300) , ENT( 300) 

COMMON  /C6TEMP/  ET,  P7 

COMMON  / B5DAT A/  V0(9) ,A1, A2 , A3, DV( 9) , T AUO, NSA ( 300 ) , PM ,GAMM1 v 9 ) , 

1  FRACT2( 300 ) ,Y1( 300) ,XE0(300) ,VP,V2 
VP  =  SPECIFIC  VOLUME  AT  WHICH  HUGONIuT  INTERSECTS  PHASE  BORY 


INTEGER  H, BURN, S, SI, ZON, CYCLE, COUNTS, CYCLES, ALP, OPT  ION, H2,HS1,HS, 
l  BURNS, HS2 
REAL  L,M, LINEAR, LEFTP 


C 


GO  TO  14 

ENTRY  B5EQST( S, J ) 
GO  TO  121 


o  n  o  o  f‘000 
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ENTRY  POINT  TO  INITIALIZE  BURN  $*****************¥****************** 


14  RHO ( S ) =7 . 84 
A 1= 1 . 667 
A2=3.4 
A 3=0.0 
DV(SI*-.004 
PM=. 130 
CVl=.45E-05 
CVMIX=.46E-05 
GAMM1 ( S ) =1 .6 
F0=0.0 
T0=300.0 
DPDTHX=-6.5E-05 
TAUO=3.0 

VO(S)=1.0/RHO(S) 

VP=VO(S)  /(!  .0*(-Al+SQRT(Al**2.*4.0*A2*PM)  )/(2.0*A2)  ) 

V2=VP*DV( S ) 

WRITE (6, 960)  RHO( S ) « A  1, A2 , A3 ,DV<  S ) , PM, CV1 , DX «  S) , CVMI X ,GAMM1 { S ) t  EO, 
1  TO, DP DT MX,  TAUO, H (S)»L(S) 

CSPS=.5 
HSl=(H(S-l)+l) 

HS2=H(  S )  *2 
DO  39  J=HS1, HS2 
V(J)=VO(S) 

33  U(J+1)=0.0 
0 ( J ) =0.0 
P ( J ) =1.01346-6 
TL IHAI J )=DELT 
CSP(J)=CSPS 
V1(J)=V0(S) 

FRACT2I J )=0.0 
E( J)=0.0 
ENT ( J ) =0 .0 
TMP ( J  )=T0 
NSA( J )=l 
XEQ(J)=0.0 
39  CONTINUE 
RETURN 

ENTRY  POINT  TO  SET  EQUATIONS  OF  STATE  FOR  RURN5 ********************* 

121  N$W=N$ A( J I 

GO  TO  (220, 22?), NSW 

C - MATERIAL  IS  IN  «  MASE  1 

220  FTAMl=(VO(S)/VN)-l.O 
PT=A1*ETAMI+A2*ETAM1**2 

CSP(J)  =  ( Al*VO(S)*2.*A2*VO( S)*( VOC  S) /VN-1* ) *3. *A3*V0( S) * C Vn( S ) / VN- 

1  l s )**2. )**. 5 

ET=E( J)-0.5*(P( J)*PT*QA*Q(J1)*(VN-V( J) ) 

IF(ABS(PT).LT. 1 .OE-51 PT=0.0 
IF(PT.GE.PM)CALLZMIX( S,J) 


o  o  o 
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PI JI=PT 
E l J )  =  ET 
RETURN 

? 22  CALL  ZKIXIS,J) 

E ( J )  =  ET 
PI  J)=PT 
RETURN 

960  FORMAT{1HO,5X,3HRHO,12X,2HA1,12X,2HA2,12X,2HA3,12X,2HDV,12X,2HPM, 

1  11X,3HCV1,12X,5HDXIS ), /I  X,  8E14.  6, /1H0, 4X, 5HC VMt X ,9X , 5HGAMMI , 1 1 X, 

2  2HE0»  12X,2HT0,10X, 6HDPDTMX, 11X,4HTAU0 , 11X,4HHI S) , 1 l X ,4HL I S ) » / 

3  6F14.6,4X,I7,3X,E14.6) 

END 

SUBROUTINE  ZMIX(S.J) 

THIS  SUBROUTINE  SUITABLE  FOR  COMPRESSION  PHASE  ONLY 

COMMON  /C1Z0N/  HI  9 ) ,  BURN I9)»L(9),DX{9) »Sl»  RHO (9 ) 

COMMON  /C2TIME/  T IMES .CYCLE, DELT, DTNtD TMX , TLI MAI  300 ), JCR! T , 

I  TQUI'.tTAU 

COMMON  /C4FL0H/  U 1 300), VI  300), XI  300 ) ,  Q I  300 ) , P 1300 ) , E I  300 ) »QA,VN, 

1  MASSI300),CSPI300) 

COMMON  /C6TEMP/  ET,  PT 

i  COMMON  /B5DATA/  V0( 9 1 , Al ,  A2 , A3 , DV l 9 ) , T AUD, NSA I  300 > , PM ,GAMMH 9 ) , 

1  FRACT2I300»,V1( 300) ,XEQ(300) ,VP,V2 

C 

INTEGER  H, BURN, S, S 1, ZON, CYCLE, COUNTS,  C YCLES, ALP, OPT  ION, H2,HSi,HS, 
1  BURNS, HS2 

REAL  L, MASS, LINEAR, LEFTP 
C 

1  C 

'  NSA( J  5=2 

XQ=FRACT2( J ) 

X  EQO=XEQ I J ) 

CA=TAUO*DELT 
IFIVN.GT.VP)  GO  TO  2 
IF  (VN.GT.V2)  GO  TO  3 
XFON= 1.0 
GO  TO  6 
2  X EQN=0 .0 
GO  TO  6 

f  3  XEQN=1.0+I VN-V2»/DVISI 

I  6  CONTINUE 

XN=(XO*l l.O-CA/2. 0) +0.5*1 XEQO+XEQN) *CA ) / 1 1 . +C A/2. 0) 

IFIXN.LT. 0.0)  XN~ 0. 0 
Vt=  VN-XN*DV C  S ) 

EMU1=(V0(S)/VT)-1.0 

i  ‘  PT=A1*EHU1+A2* EMU  1**2. 

CSPIJ)=I  Al*VOIS)+2.*A2*VOIS)*( VOI S) /VN-1. I +3. *A3*V0t SI  * IVuI S) /VN- 

f  l  1  * )**2. )**. 5 

E  ET=EI  J)-0.5*(PT+P(J )  +  QA+Q< J) J  +  IVN-VI J) ) 

(  V 1 I J )=VT 

\  FRACT2IJ)=»XN 

t  XEQI J )=XEON 

RETURN 
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END  ' 

SUBROUTINE  WRITE 

COMMON  /C1Z0N/  H< 9) , BURN? 9) , L ( 9 ) ♦ DX( 9) , SI ,RH0(9 > 

COMMON  /C2TIME/  T IHES  tCYCLE,  DELT,  DTN,DTMX,  TL I  MA(3Q0)  ,  JCRI T , 

1  TOU 1 Tf  TAU 

COMMON  /C3CTRL/  COUNTS, JS TAR , JPE , JPB, JQU1 T, LAST .CYCLES 

COMMON  /C4FL0W/  U<  300 1 ,  V(  300 1 ,  X<300),  QI  300 ) ,  P  (300  ) ,  E  (  300  J  ,QA,VN, 

1  MASS(3001,CSP(300) 

COMMON  /C5THFR/  TMP  <  300) , ENT  I  300) 

COMMON  /C7GNRL /  ALP .OPTION, CONA, CQ, LEF TP 

COMMON  / B5DAT A/  V0( 9 5  , A1 , A2 , A3, DV ( 9 ) , TAUO, NSA ( 300 ) , PM , GAMM1 ( 9 ) , 

1  FRACT2 ( 300) ,V1( 300) , XEO ( 300 ) , VP , V2 
C 

INTEGER  H,BURN,S,S1, ?0N, CYCLE, COUNTS, CYCLES, ALP, OPT  ION, H2,HS1,HS, 
1  BURNS, HS2 
C 

REAL  L, MASS, LINEAR, LEFTP 
C 

GO  TO  14 


ENTRY  WRITE1 
GO  TO  121 
14  WRITEI6,  302  ) 

WRITE (6,  304 ) T IMES, DEL T,DTN, CYCLE, JCRI T 
WRITE (6, 306) 

S  =  2 

I  F  (  JPB,  EQ.  1  .AND.  JSTAR.  GT.H  1 2 )  *-10  )  JPB=H(2) 

DO  330  J=JPB, JPE 
IF(J.GT.HIS) )  S=S+1 

?10  WRITE  I  6, 318) J,UIJ),V( J),P(J),E(  J) ,QIJ) ,FRAC T2 < J ) ,  VI  (  J ) , X(  J) , TMP { J ) 
l , TLIMAI  J ) 

330  CONTINUE 

C - NEXT  TWO  STATEMENTS  (COMMENTS)  ARE  TO  BE  USED  IF  GRAPHING  IS  DEStRED 


CALL  MANUAL! 1.25*LEFTP,0.  ) 

CALL  GRAPHKP, JPE) 

IF(LAST.EC.1)CALL  EXIT 
RETURN 

121  WRI TE( 6, 862 ) 

DO  46  J  =1,2 

46  WRITE(6,962)J,U(J ) , V( J ) , P ( J ) , E< J ) , Q! J )  ,  FRAC T2 ( J )  ,  VI  (  J ) ,X(  J)  ,  T 

-MP( J ) ,TLIMA(  J ) 

DO  57  S=2» S 1 
HS1  =  H( S i —  1 


HS2“H( S ) +2 
DO  57  J=HS1,  HS2 

57  WRI  IE  ( 6,  962  )J,U(J',V(J!,P(J),E(J)  *  Of  J)  *  FP.ACT2 1 J  5  r  VI !  J )  ?  X (  J)  •  T 

-MP( J  )  ,TL IMA( J ) 

RETURN 

302  FORMAT! 1HI ) 

304  FORMAT (10X,6HTIME=  , E  14.8 ,4X, 6HDEL T=  , E14.8,4X, 5HDTN=  ,E14.8,4X,7H 
-C YCLE=  , I5,4X,7HJCRIT=  ,I5/5X) 

306  F0RMAT(2X,1HJ,6X, IHU, 9X, 1HV, 9X, IHP , 9X, 1HE ,9X , 1HQ , 7X , 6HFRACT2 , 6X , 

1  2HV1»9X, IHX, 9X , 3HTMP, 6X, 5HTLI MA//5X ) 

318  FORMAT(I4,8F10.6,F7.1,E13.5) 

462  FQRMAT(2X,1HJ»6X, IHU, 9X, 1HV , 9X , IHP, 9X, 1HE , 9X , IHQ, 7X , 6HFR ACT2 , 6X , 

1  2HVI,9X,1HX,9X, 3HTMP, 6X, 5HTL I HA//2X ) 
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962  F0RMAT(I4,8F10.6,F7.1,E13.5) 

3  F0RMAT(1H0,2X,6HTIMF=  , E l 4. 8 , 4X , 5HUFS=  ,  F 14. 8 ,4 X, 5HXF S  =  ,E14.8,4X, 
-7HCYCL  E=  , I4,3X,7HJPMAX=  .  I  4, 3  X, 7HJCRI T  =  *  1 4/ 2X  > 

END 

SUBROUTINE  GRAPH4I A ,B  ,C,D ,N> 

COMMON  / C3CTRL /  COUNT S, JSTAR , JPE , JPB, JQUI T, L A  ST ,LYCL ES 

COMMON  /C4FL0W/  U  (  300  J , Vt  300  ) , Xt 3001,3(300) , P (300) ,E ( 300) ,OA,VN, 

1  MASSC300)  ,CSP(  300) 

COMMON  /C7GNRL/  ALP, OPTION, CONA,CQ,LEFTP 
REAL  L, MASS, LINEAR, LEFTP 

OIMFNSION  A(N),  B(N),  C(N),  D(N),  P0INT(4),  GRAPH(122) 

DATA  POINT/IMA, 1HB, 1HC, 1H0/ 

OATA  MSWTCH, BLANK, PERI0D/0,1H  ,1H. / 

M=4 

I F ( MSWTC  H« EQ . 1 )  GO  TO  50 
AMAX=D ( 1 ) 

AMIN=D( 1 ) 

00  1  1  =  1, N 

IFIAMAX.LT.Ot  I ) )  AMAX=0(t  ) 

1  I F( AMIN. GT.D( I ) )  AMIN=D(I) 

GO  TO  2 

ENTRY  GRAPH3(A,B,C,N) 

M  =  3 

IF(MSWTCH.EQ.l)  GO  TO  50 
AMAX=C ( 1 ) 

AM IN=C ( 1 ) 

2  DO  3  1=1, N 

I F ( AM AX.LT .C (  I  )  )  AMAX=C( I  I 

3  I F ( AMIN. GT.C( I ) )  AMJN=C(I) 

GO  TO  4 

ENTRY  GRAPH2( A,B,N) 

M  =  2 

IFIMSWTCH.EQ.'l )  GO  TO  50 
AMAX=B ( 1 ) 

AMI N=  B ( 1 ) 

4  DO  5  1=1, N 

I F ( AMAX.LT  .8(1))  AMAX=B(I) 

5  I F( AMIN. GT.B( I ) )  AM  I N=B< I ) 

GO  TO  6 

ENTRY  GR APH1 ( A, N ) 

M  =  i 

I F ( MSWTCH. EQ , l )  GO  TO  50 
AMAX  =  A( 1 ) 

AMINsAt 1) 

4  DO  7  1=1, N 

IF ( AMAX.LT. A ( I )  »  AMAX=A( I ) 

7  I F ( AMIN. GT. A( I ) )  AMIN=A(I) 

50  SC=( AMAX-AMIN1/120. 

WRITE (6, 100)  AM IN»AHAX, SC, (PERIOD, 1=1, 121) 

100  FORMAT ( 1H1, 17HRANGE  OF  GRAPH  IS,F15.8,8H  THROUGH, F15.3//1X*17HSCAL 
11NG  FACTOR  IS, F15 .8//  105X,21H1U  111111  Ulllllll 111 /15X.11 1H111 1  111 

21 1122222222223333333333444444*444555555555566666666667777777777088 
38  88888899999999990000000000111 llill 112 /6X,120H12345678901 23  V, 67890 
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4123456789012345678901234567890123456789012345678901234567890123456 
578901 2345678901 2345 678901 234567890/5X,  121 A1 1 
00  8  I a  1 ,  N 
DO  17  Ja 1, 121 
17  GRAPH ( J 1 =BLANK 

GO  TO  (9,10,11, 121, M 
1?  K4-(D(l»-AMIN)/SC+l. 

I  FU4.LT.0.0R.K4.GT.122)  K4=122 
11  K3»(Cm-AMIN)/SCn. 

IFIK3.LT.O.OR.K3.GT.122)  K3=122 
10  K2=(B(I)-AMIN)/SC+i. 

IF(K2.LT.0.0R.K2.GT.122)  K2=122 
«  K 1= « A(l)-AMIN)/SC*l. 

I F ( K 1 ,LT .O.OR.KI.GT.122)  Kl= 122 
GO  TO  (13, 14, 15, 16), M 
16  GRAPH(K4)=P0INT(4) 

15  GRAPH«K3)  =  P01NT(  3) 

14  GRAPH (K2 )  =  P0 INT ( 2  ) 

13  GRAPH (Kl)=POINT(ll 

IF(GRAPH( 121) .NE. BLANK)  GRAPH (120) -GRAPH ( 121 ) 

WRITE (6, 101)  I, l GRAPH! Ill ,11*1,120) 

101  F0RMATI2X, 13, 1H.,  120A1) 

8  CONTINUE 
MSW1 CH=0 
RETURN 

ENTRY  MANUAL ( A 1 , A 2 > 

AMAX=A 1 
AM IN*A2 
MSWTCH=l 
RETURN 
END 
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For  numerical  output  of  this  problem,  refer  to  Vol.  II  of 
this  report  filed  in  the  Document  Library  of  Ballistics 
Research  Laboratories. 
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D.  List  of  Labels 
DECIDE 

DX(S)  =  Eulerian  space  interval  in  region  S  at  t  =  0 
=  L(S)/Z{*N(S) 

H(S)  =  no.  of  cells  from  left  boundary  through  region  S 
S 

=  y  Z($N(L) 

L~> 

L=2 

B5INIT(S) 

RH0(S)  =  density  at  zero  pressure  in  region  S 
Al,A2,A3  =  coefficients  in  Eq.  (4.5) 

DV(S)  =  v2(p,T)  -  v1(p,T) 

PM  =  pressure  at  which  the  Hugoniot  in  phase  I  intercepts 
the  phase  boundary 
CVl  =  Cvl 
CVMIX  =  CVIn 
GAMMl(S)  ■=  r 

E0  =  internal  energy  at  the  foot  of  the  Hugoniot 

Tf)  =  TQ 

DPDTMX  =  (^p/^T) v 

TAU0  =  1/t,  Eq.  (5.11) 

VP  =  specific  volume  in  phase  I  at  p  =  PM  =  v^(pM,T) 

V2  «  v9(pM,T) 

CSPS  =  starting  value  for  sound  speed 
J  =  index  for  space  grid 
V(J)  =  vj 

U(J)  =  U. 

.1 

Q(J)  “ 


I 
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p(j)  -  Pj 

TLIMA(J)  =  value  of  Atj  for'  next  time  step 

CSP(J)  =  sound  speed  in  cell  J 

E(J)  =  Ej 

ENT(J)  =  Sj 

TMP(J)  -  Tj 

NSA(J)  =  switching  index 

=  1,  phase  I 

=  2,  mixed  phase 

=  3,  phase  II 

MAIN 

X(J)  =  Xj  (Fig.  3.2) 

MASS(J)  =  mass  of  cell  J 

JSTAR  =  cell  label  just  ahead  of  shock  front  at  which 
computation  stops  for  each  time  cycle 
TIMES  =  t 

CYCLE  =  number  of  times  t  has  been  incremented 
JCRIT  =  value  of  J  for  which  TLIMA  is  minimum 

LAST  =  switching  index  for  halting  program  after 

writing  last  output. 

PPEAK  =  maximum  computed  pressure  in  each  cycle 
TLIMB  -  TLIMA  (JCRIT) 

PLEFT  =»  pressure  applied  to  left  boundary 
DFNU  =  mass  in  cell  J+l 
XA  =  x(t  +  At) 

VN  =  v(t  +  At) 

QA  =  Q(t  +  At) 

JPMAX  =  value  of  J  at  which  p  is  maximum 


ZMIX 


FRACT(J)  =  cvj 

XE‘4(J)  = 

Vl(J)  =  vXj(p 


i 
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Fig.  1 .  --FLOW  CHART  FOR  BURN 
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r 


Q  U(JSTAR  +  1)  >  10"5  n) 


JSTAR  = 

FSTAR  +  1 

- x 

* 


CYCLE  =  CYCLES 
or  TIMES  &  TQUIT 
or  J  =  JQUIT 
© 


MOD  (CYCLE,  COUNTS)  =  0  |y} 
- — - 


40 


LAST 


WRITE 

PRINT  FLOW  VARIABLES 


n  LAST  =  1  y 


Return 

to 

MAIN 


CALL 

EXIT 


I 


Reset  DELT 
Go  to  40  for 
new  cycle 


Fig.  1.  (c) 
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A  procedure  is  described  for  developing  simple  approximate  equations  of  state 
of  liquids  from  Hugoniot  P-V  relations  determined  In  shock  wave  measurements. 

This  Is  applied  to  a  number  of  liquids  and  a  table  of  coefficients  is  given. 

The  formalism  of  irreversible  thermodynamics  is  applied  to  time-dependent 
phase  transitions  In  Iron  and  an  approximate  set  of  constitutive  relations  Is 
obtained  In  a  form  suitable  for  numerical  integration  with  the  equations  of 
continuum  dynamics.  These  are  applied  in  an  approximate  form  to  study  the 
development  of  the  two-wave  structure  in  Iron  caused  by  the  <*-e  phj.se  transition. 

Finite  strain  theory  Is  applied  to  the  analysis  of  shock  wave  data  for  quartz, 
and  the  results  supply  enough  information  to  estimate  some  of  the  fourth-order 
elastic  constants. 
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